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SUMMARY SYMBOLS 


An equation for estimating wing weights is developed which h = wing span, ft 
gives reasonably accurate results Fig. 1 shows a comparison Cr = root chord, in 
of predicted weights with actual weights of a large number Cy = tip chord, in 
of aircraft varying from supersonic single seat fighter types to dp = maximum thickness at root chord, in 





multiengined low subsonic speed aircraft. Military as well as dy = maximum thickness at tip chord, in 
commercial types by eleven major manufacturers are included f = equivalent stress, Ib. per sq.in 
Many more aircraft than those shown on Fig. 1 were examined a stress analysis gross weight, Ib 
but not included because they were earlier or less developed ver h = fraction of semispan measured from tip 
sions of the particular ones that are shown n = ultimate design maneuver load factor 
The variables included in the equation consist of the wing P= concentrated load, Ib 
geometry dimensions, basic flight design weight, ultimate load = Colt 
factor, maximum equivalent sea level speed and an equivalent = maximum equivalent level flight speed at S.L., knot 
stress factor, f. The values of f may be determined from Fig. 3, ‘. = total wing weight, Ib 
together with the recommendations given on page 372 wp = air loading at root, Ib. per unit area 
The basic wing weight equation is wy = air loading at tip, lb. per unit area 
Ww = (GQ: — Q2) Wp, = inertia loading at root, lb. per unit area 
Wy, = inertia loading at tip, lb. per unit area 
where ; 
a (dp — dr)/dy 
14.4b6n 3] B = thickness ratio 


¢ angle of sweepback of 1/4 chord In 


T ¥ . - . . 
density of material, lb. per cu.in 


14.460n 


O.O0lrea + O.257r 0.0958 + 0 105 


2.4nb ‘ ( 3b j ) PC] s/i 
= 7 ae ( — ‘) = 
ad7 cos ¢ 7 
’ 12h + + 0.0958 + 0.105 
P ( + i) y; 
cos ¢ ‘ " 2) 
7 P Pp 


and the remaining symbols are defined in the next section. Al . , 
though this equation appears to be complex, a wing weight can be 12h 
calculated in approximately 30 min. The weight will include >» P( 
flaps, ailerons, as well as leading-edge slats or flaps in the case of od ne 
swept wings. The effects of distributed wing weight and concen att 2. ob te £4 oe Bs (2 +7) + 21 + 2 
trated weight items in the wing are taken into account distance of spanwise center of pressure from the plane 
of symmetry of the wings 

= (tan Va — tan! Vha)/Va 

V a 

1] 
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Subscripts | RE | eh 
d = distributed FIGU 
piped ssiccsesile COMPARISON OF ACTUAL VS COMPUTED 
R = root WING WEIGHTS 
T = tip 30 
=o 
INTRODUCTION _j20 
A OVERESTIMATE OF THE wing weight of anairplane 615 
in the preliminary design stages (or too high a ie) 
weight estimate of any other part of an airplane) will + 10 
cause the resulting airplane design weight to be too E 8 
high by a factor of perhaps 5 to 15 times the amount © 6 
by which the wing weight is too high. Since the wing = 5 
weight accounts for approximately 10 per cent of the o 4 
airframe weight, it can be readily appreciated that ac- Z 
curate wing weight estimation is indeed a very impor- = 3 
tant element in airplane design. | 
The design of an airplane consists of a process of DD 2 
iteration. Before any airplane design is selected for o 
detailed engineering analysis, a number of preliminary © 
designs are made and the effects of varying airplane | 
: ts igs: | 2 3456 81 °5 20253 


parameters are investigated with the purpose of ob- 
taining the best combination. To be generally useful, 
therefore, a method of wing weight estimation must be 
accurate, rapid, and be based on such airplane dimen- 
sions and parameters as are easily determinable in the 
preliminary design stage. 

Because most wing weight estimating methods have 
empirical constants based on data from existing air- 
planes, work on improved methods reflecting informa- 
tion from the more recent airplanes is continuously in 
progress. This report was undertaken to revise and 
improve an earlier wing weight-estimating method 


prepared by the author. * 


ANALYSIS AND DERIVATION 


A modern airplane wing is a cantilevered structure 
which is subjected to shear, bending moments, and tor- 
sional moments resulting from the applied air loads. 
The wing strength must be adequate to withstand the 
applied forces. Therefore, it is natural to relate the 
wing weight to the strength required, provided the 
method is sufficiently simple and accurate. 

The examination of a large number of wings has re- 
vealed that, in general, the strength required and, con- 
sequently, the weight, increases rapidly from the tip in- 
ward. This seems to be true, not only for the basic 
wing structure, but also for the ribs, flaps, ailerons, and 
slats. In other words, a plot of cross-sectional wing 
weight per inch of span versus spanwise location would 
show a rapid rise from the tip inward. 

The above analysis led to the hypothesis that the en- 
tire wing cross-sectional weight, that is, not only main 


* NavAer Report No. DR-1139, “Aireraft Design Analysis 
Method as Employed by the Research Division of the Bureau of 
Aeronautics, U.S. Navv Department,’ by Ivan H. Driggs, Octo 
ber 6, 1949; Journal of the Royal Aeronautical Society, Vol. 54, 


February, 1950. 








COMPUTED WING WEIGHTS - OOO LB. 


structural items, but secondary structure and other 
items included might be proportional to the shear, and 
bending moments, since for a cantilever beam, the shear 
and bending moments increase from the tip to the root. 
The torsional moments were eliminated from consider- 
ation because of the complexity of the resulting equa- 
tions and the unavailability of the necessary torsional 
loading parameters in the preliminary design stage. 
The main assumption in the derivation of the wing 
weight equation is, therefore, as follows: 

(1) The weight per unit span of any chordwise cross 
section of a wing is proportional to the shear and bend- 
ing moment forces at that cross section. 

Additional simplifying assumptions used are: 

(2) The air loads and wing structural weight vary 
linearly from tip to root or are uniformly distributed 
over the entire wing area. 

(3) The wing is assumed to be a box with centerline 
and tip cross sections being rectangular and having the 
long dimension equal to the respective chords and the 
short dimensions or depth equal to the respective maxi- 
mum thicknesses. <A linear variation of chord length 
and thickness is assumed between € and tip. 

Other less important assumptions will be stated in 
the derivation. The manner in which torsional effects 
are brought to bear will be discussed later. 


Derivations 





Let 

a= (dp = dr) dr; r= Cr Cr and w = Wr/ W7 
We and wr are the air loads in pound per _ unit 
area, 
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A METHOD FOR 


The spanwise center of pressure is given by 
¢= [w(1 +7) + (1 + 37) + 
[2w(2 + r) + 2(1 + 2r)] (1) 
zis measured from the € of the wing. 
From Case I, Diagram I: 


= A + Bx (2) 


Wr = wr + (2x/b) (we — wr) 


C, = Cr + (2x/b) (Ce — Cr) = C+ Dx (3) 
d, dr + (2x/b) (dr — dr) = dr{1 + (2ax/b)] (4) 
The shear at station x can be determined to be 
F, = (x/6) (2QwrCr + w,Cr + wrC, + 2w,C,) 
= (x/6) [(6AC + 3(BC + AD)x + 2BDx?] (5) 


and the bending moinent at x, 


12) (3wrCr + w,Cr + wrC, + w,C,) 
12) [(6AC + 2(BC + AD) x + BDx?] (6) 


M, = (x? 


= (x* 


If it is assumed that the bending material is concen- 
trated at the upper and lower surfaces of the box of 
Diagram I,* and if the compressive and tensile stresses 
are equal, then the area of bending material at any sta- 


tion of one side of a wing is given by 


A a= 2M,, (d,fx) (7) 


and the area of shear material is 
As =P.) fs (S) 


A limited investigation was made of the mean stress 
of all material put into a wing to resist shear stress pri- 
marily and compared to the mean stress of all the ma- 
terial in a wing carrying bending stress primarily. It 
was found that the bending stresses were of the order 
of 4 times the shear stress. The assumption is there- 
fore made in this report that 


* Slope of upper and lower surfaces is neglected. 
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S 
ar DIAGRAM I 
CASE I- LINEARLY VARYING LOADING 
fp = tf 
Hence, from Eq. (4) 
A s=4 F,, ‘te (9 ) 


Another investigation of a number of aircraft was 
made to determine the variation of bending stress with 
span, assuming that all the bending material was con- 
centrated at the upper and lower surfaces of a box, as 
shown in Diagram I determined by the wing geometry 
of the aircraft chosen. The results are shown in Fig. 
2. Also shown in Fig. 2 is the parabolic spanwise vari- 


ation assumed in this analysis, that is, 


Vax = fV2x/b (10 


fry - 


where f is the equivalent stress at the € of the wing. 


With the assumptions made thus far, the total wing weight, due to distributed air loads, is given by 


"hb 2 b/2 
Ww, = 2 | pAm dx + 2 I pAs dx 
0 0 
b/2 DAF, 
dx + 2 —— dx 
0 te 


II 


*b/2 9 
9 / p2 ate 
0 df 


(11) 


Substituting Eqs. (3), (4), (5), and (9) into Eq. (11) gives: 


1 , f - 2p(6A C+ 2(BC + AD)x + BDx?] (<*, 12) F 
a = 2 sl 
0 


dr[1 + (2ax/b)] f V2x/b 


p(b/2)' 
3f 


1 + (2ax/b) 


dr 


> ] 


ie 


' f > 4p[6AC + 3(BC + AD)x + 2BDx’] (x/6) F 
2 — ax 
0 | V 2x b 


i] l f 2 (BDx'* + 2(BC + AD)x** + 6ACx”?] 
——" dx +- 
0 


> 


6/2 
4 f (2BDx*? + 3(BC + AD)x** + 6ACx'*] ax | (12) 
0 
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Eq. (12) may be solved to give 


’ 2pb*Crwr jb [(w—-— 1) 1-7) .. (w—-1)+(1—r) .. : 
Ws. * can ‘| ; Y, + : y, + | + 
0.105 + 0.095 + 0.257r + 0.017, (13 
where 
Yo = (1/7a@) — (1/5a?) + (1/3a%) — (1/a*) + [(tan-! V a)/a *] 
Vi = (1/5a@) — (1/38a?) + (1/a*) — [(tan-! V a) a’*] 
Y, = (1/8a) — (1/a?) + [(tan-! V a)/a’’] 


It will now be assumed that the wing dead weight is distributed in the same manner as the air loads or 


(Wr WT )air loads = (Wr, WT 1) dead weights — 


9 


Hence, the wing weight reduction due to distributed wing weight is given by an equation analogous to Eq. (13 
with wy replaced by w7,. 
The weight due to dead weight or inertia loads would be 


Way = [(2pb?Crwr,)/f] Z 


where Z is the quantity in {| | of Eq. (13). 
The wing weight due to distributed loads becomes 


Wop = Wea ing ioe 
= (wr — wWr,) [(2pb*Cr Z)/f]. 14) 








FIGURE 2 
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FIGURE 3 
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‘yy. /100 
v, /' 


The total air load on the semiwing is given by 

nG = (b/6) (2wrCr + wWrCr —- wrCr = 2wrCr (15 
and the total inertia load can be written as 
nW,* = (5/6) (2wr,Cr + Wr C7 — wr),Cr + 2wr Cr) (16 
Subtracting Eq. (16) from Eq. (15) gives 


n(G — W.) = (Wr — WT) [(b 6)Cr] [(2 + Q)r + (1 + 2w) | (Wy — Wr,) = 
n(G — JW?) /[(b/6)Cr] [(2 + w)r + (1 + 2w)] (17 


By eliminating (wr — wy,) from Eq. (14), the wing weight due to distributed loads, 
3 12pbn(G — W,,) bif@—1l)al—-r).. (w—1)+ (1 — 7) : 
W., = . Yo + Y,+ r¥2 <4 
tHrw + 2r + 2w + 1) l4dr 6 BS 
0.01lrw + 0.257r + 0.095wH + 0.105 (1S 


Eq. (18) is in consistent units. Normally } is given in feet whereas dr is usually given in inches. The density 
of aluminum p = 0.1 Ib. percu.in. Also, in the case of angles of sweep of the 1/4 chord line greater than 15° the 
span will be divided by the cosine of the angle of sweep (an assumption based on the fact that the true length of 
the wing is of most interest rather than the normal span). With the incorporation of these items, Eq. (18) becomes 





* The entire wing weight is used, that is, the effects of concentrated loads are included, since the source of the structural weight 


is immaterial 
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; 14.46n(G — W,,) j 3b 
Wop = 


» Vy + 
f(rw + 2r + 2w + 1) cos ¢ (dr cos ¢ 


6 3 


(w — 1) (1 — 7) ‘ (w — 1) + (1 — 7) ; 
‘ y,+ ry, + 


0.01rw + 0.257r + 0.095w + 0.105} 19 


Eq. (19) gives the wing weight when the wing does not carry any relieving concentrated loads, that is, |’ » be 
comes equal to IV’,. 


When 
~— 14.4bn j 3d E —1)(1—-~+r) y, + (w—1)+ (1 —~v7r) y+ ry] 4 
2 (rw + 2r + 2w + 1) cos ¢ ldr cos ¢ 6 3 
0.01rew + 0.257r + 0.0950 + 0.105, 
Wun = [((G — W,) Q1]/f 20 
Now when IV’, = IW, and solving Eq. (20) for IV, 
We = (GQ)/U + Qi) (for distributed loads only) 2 


The effect of concentrated loads will now be considered. From Diagram II (opposite page), 


d, = dr + [(2x/b) (de — dr)] = dr} 1 + [(ax) b]} (22 
M, = —nP|x — h(b/2)] (23 
F, = —nP 4 


Then the area of material required at any station to withstand bending is 


| 2nP|x — h (b/2)] 
d, les 
and shear area 


A, = —[(4nP)/fs,| lsee Eq. (9)] 


The wing weight due to a concentrated load is 


> 


; “62 2QnP[x — h(b/2)] “6/2 AnP 
W.. = —2p dx — 2p dx (25) 
JS hb 2 d,fr, h(b’2) Ser 


Combining Eqs. (10), (22), and (25) 


WV tpn P(b/2) | l I. 2 [x * — h(b/2)x |] dx 42 ba ] - 
a é x x (26) 
/ dr Snio 2 }1 + [(2ax)/b]} nine ~ 


let 
y? = (2anx)/bd; 2y dy = (2a/b) dx 

x’? = (b/2a)"y; when x = 0/2, y= Va; when x = h(b/2),y = Vha 

iW, | b? F i « y? dy hb [- « dy ‘ 2b | i ‘ ™ 
= — = ay (24 

tpnP 2dr a V hea l + Se 2dra : ~“vVha I a ¥* a Vha i 
. tpnPb tan-! Va — tan-! Vha 

W.. = - | (6 + 4adr) (1 — Vh) — (6 +h) - (28 

2ad7 Va 


when 4 is in feet and p = 0.1 Ib. per cu.in. and dividing 6 by cos ¢ 


2.4 nPl 3 121 tan-? Va — tan Vha\* 
= = = | ( + adr (l-— Vh) - ( + n) = = — = (29 


adyf cos ¢ cos ¢ 0S © 
7 ¥ L ¥ C ¥ Va 
When concentrated loads occur at more than one station along the span, then 


* When a is negative, the quantity within the parentheses is determined without regard to the negative sign of a 
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v.«- rie : 4 (= + air) > Pa — Vh) — SP (— + h) vs | (30) 
where 
Y, = (tan-!' Va — tan! V he) Vea 
When 4 = 0. 
Ws. = —[(9.6nbP)/(f cos ¢)]} [3b/(dr cos ¢)]V¥3 + 1} (31) 
Let 
hm = me : E See + adr) oP - Va) - XP (= re h) | 
The wing weight of a wing with distributed and concentrated loads will be 
We = Wen + Wee = {1G — W.)Q)/f} — (Q2/f) 
= (GQ; — Qe)/(f + 1) (32) 


Case II is the same as Case I except Cr = 0 and dr > 
0. Eg. (32) is directly applicable by setting Cr = 0 
wherever it occurs. 

Case III is the same as Case I except dr > dr. Eq. 
32) is directly applicable. The values of a will be 
negative but }4, Yo, and V; will be positive and can be 
determined from Fig. 4. 

Case IV is the same as Case I except Cr > Cr. Eq. 
32) applies with proper regards for signs. 

Case V is the same as Case I except Cr = Cr. Eq. 
32) applies. 

Case VI is the same as Case I except for the uniformly 
distributed air load. Eq. (32) applies by setting w = 1. 

All the terms of Eq. (32) are available directly dur- 
ing a preliminary design stage of an airplane except the 
terms, }y, V1, Yo, V3, Ys, and f. The values of Yo, Vi, 
Vo, ¥3, and }y have been computed for a broad range of 
values of a, which in turn is a function of root and tip 
chords and thickness ratios, and are plotted in Figs. 4 
and 5. The values of f can be determined from Fig. 3. 
The discussion part of this report should prove helpful 
in adding judgment to the choice of the most applicable 


value of f. 


DISCUSSION 


In order to determine the value of /, the character- 
istics of more than 50 aircraft were obtained and Eq. 
Unfortunately, the spanwise 


32) was solved for f. 
center of pressure of most of the aircraft was unknown 
and consequently a uniform loading was used except in 
some 5 or 6 cases. The values of f were then plotted 
against (wG)/\y. Inspection of the resulting graph 
indicated that although the spread of points was large, 
the latest models of different designs and different 
manufacturers showed improvement over earlier models 
and appeared to be grouped closely about a mean 
line. 


The various aircraft were then grouped according to 
manufacturer and type and plotted as shown in Figs. 
6 and 7. It became apparent at once that the value 
of f, in most cases, appeared to follow some learning 
law. In the majority of the cases the value of f im- 
proved markedly for the production model over the ex- 
perimental model. The figures show ten cases (A, 
B, D, E, F, G, I, K, N, and U) of series of airplane 
models beyond the first production model. In nine 
of the ten cases, succeeding models showed improve- 
It is suspected that the situ- 
Very likely the 
data that were available for analysis are incorrect. 

The U series of aircraft tend to bear out this con 
An analysis of the series U from data that are 


ment over prior models. 
ation is true in the tenth case, E, also. 


tention. 
normally available in stress reports and weight reports 
gave the results shown on Fig. 7. Later in the investi- 
gation very detailed analysis became available about the 
composition of the wing weight as given in the official 
weight report. The data showed that several items of 
weight for provision of landing gear cut out, wing fold, 
tip tanks, etc., could be properly eliminated and a net 
structural weight determined. Values of f were recal- 
culated with the new weights and the results are shown 
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is the series U’. The reversal of trend as shown in the 
U series of airplane model 3 has been eliminated and the 
None of 


series U and U’ has been plotted on Fig. 3 because the 


value of f increases steadily. the models of 
inboard region of the wing deviated too greatly from 
the values used in the calculations. The analysis, how- 
ever, does properly reflect the change of trend between 
the U and U’ series. 

The evidence of Figs. 6 and 7, representing as it does 
more than 20 series of aircraft of many aircraft manu- 
jacturers over a span of about 20 years, appears to show 
fairly conclusively that the wing structure is used more 
effectively with succeeding models of the same basic 
design. This is a natural and probably expected trend. 
During the design of an experimental model the engi- 
neer must rely mainly on analytical methods, with their 
attendant inexactness and consequent tendency to 
“play it safe,’’ to determine the strength of the wing 
structure. benefit from the 
mistakes on the experimental model and in most in- 


Subsequent redesigns 
stances are also benefited by the results of proof tests or 
destruction tests of the experimental model. 

The question still remains—what values of f should 
be used in predicting wing weights with the equations 
derived in this report? Every airplane designer knows 
that overestimating the weight of any component in the 
initial design stages of an airplane means that the air- 
plane weight is increased by several times the amount 
of the overestimate. It is also an unpleasant fact to 
stress analysts that by the judicious addition after 
static tests of, say, 5 per cent to 10 per cent of the initial 
wing weight, the wing strength can be increased perhaps 
\0 per cent. Consequently, it is the opinion of the 
author that wing weight estimates should reflect the 
upper portion of the state of the art. For this reason, 
only the models having the highest values of f of each 
series, noted by arrows in Figs. 6 and 7, were plotted 
in Fig. 3. 

The / values were plotted against the parameter (7G + 
ly. It was assumed that speed is reflected in addi- 
tional torsional rigidity requirement of a wing structure 
and sometimes the high speed gust load factor exceeds 
the maneuver load factor; hence, the greater the speed, 
the heavier the wing. 

It will be noted that three different point symbols 
This was done for the following 


reasons: a careful inspection of the wings comprising 


were used in Fig. 3. 
the group to be plotted showed that most of them had 
in addition to ailerons, flaps, or leading edge flaps or 
slats, cutouts or provisions for either one or more of the 
following—landing gear, wing fold, or engine nacelles 
which cut into the wing as contrasted to engine nacelles 
on pylons. For this group of airplanes circular point 
symbols were used. Two aircraft which differ from the 
others in the region of (1V)/ Vz below 600 in that they 
have landing gear cutouts but no wing fold are shown 


with square symbols. A uniform load distribution was 
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FIGURE 6 
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AIRCRAFT SERIES 


used in determining / for the aircraft represented by 
circular and square symbols. 

Several aircraft, O2, Q2, and V1, had exceptionally 
structures—no landing 


In addition, data 


clear, straight-through wing 
gear, wing fold or nacelle cutouts. 
for two more aircraft, H2 and W1, were available which 
indicated the penalty for either landing gear or wing 
fold in the wing weight. Furthermore, all five aircraft 
are believed to reflect the most up-to-date refinements 
in structural analysis, manufacturing techniques (ta 
pered skins, extrusions, etc.), materials, and large num 
bers of engineering man-hours. Data were also avail- 
able on the spanwise center of pressure locations of 
these aircraft. These factors are believed to combine 
to give values of f higher than the other aircraft ex- 
amined and probably the best values which, on the 
average, may be expected from the state of the art at 
this time. Eq. (32) was used to determine the values 
of f and x symbols used to plot them on Fig. 3. 

A calculation was made of the mean deviation of 
the circular points in Fig. 3 from the mean line through 
them, and it was found that the mean deviation was 
1.7 percent. The mean deviation for the crossed points 


is 5.5 per cent. 


EJGURE 7 mm PRODUCTION MODELS 


WINGS PLOTTED IN FIGURE 5 EXPERIMENTAL MODELS 
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Recommendations for Choosing Values of f 


Considering the results of the investigation as pre- 
sented and numerous calculations carried out but not 
shown, the following manner of selecting the most ap- 
propriate value of f from Fig. 3 is reeommended: 

(1) For high-speed structurally ‘“‘clean’’ wings—no 
landing gear cutout, no wing fold provision, no provision 
of equipment of any kind but including flaps, ailerons 
and/or slats—use the upper curve to determine the 
value of /. 
CP location should be reflected in the value of w in Eq. 
(32). (It is also assumed that the most diligent effort 
will be made to refine the structure as the detail engi- 


A reasonably good estimate of the spanwise 


neering progresses in order to meet the predicted 
weight. ) 

(2) For high- or low-speed wings with cutouts or 
structural discontinuities for various for 
landing gear, armament installation, cutouts for nacelles 


use the lower curve for 


purposes 


or tail pipes, and wing fold 
values of f and a uniform load distribution (w = 1). 

(3) For structurally delta 
limited amount of data indicates the upper curve for 
values of f should be used together with an estimated 


“clean”’ wings a very 


spanwise CP location. 


(4) Only major items of concentrated loads need be 


FICAL SCrirsegneses JUNE, 1954 
ment, large stores, etc., but not control systems, or 


other small items of equipment. 


Example: 
Assumed: 


G = 140,000 Ibs. 

n = 3.8 

S = 1900 sq.ft. 

b = OOH. 

Cr = 342.5 in. 

Cr = 114im.; ry = 114/342.5 = 0.333 

Be = 0.11; Br = 0.08; dr = 0.08 X 114 = 91] 
g = 40°; cos ¢ = 0.766; ccs? g = 0.588 

P, = 425 Ibs. at 50 ft. from € or h = 0 

P, = 12,300 Ibs. at 8.3 ft. from € or h = 0.834 
V~ = 600 knots. 


Also assumed is that the wing is structurally ‘‘clean.” 


Determine: a = [(dr — dr)/dr| = [(342.5 X 0.11 - 
114 X 0.08) ]/(114 X 0.08) = 3.14 

From Fig. 4, Yo = 0.0325; VY; = 0.043; Y.2 = 0.0656 

Determination of Qi: It will be assumed that the 


spanwise CP for the design condition is at 0.45 of the 
semispan. Hence, from Eq. (1): 
2(1 + 2 X 0.333)] 


0.45 = [w(1 + 0.333) + (1 + 1)] 


taken into account-——landing gear, fuel in wings, arma- w = 0.653 
, 14.4 X 3.8 X 100 ); 3X 100 —_ — 1) (1 — 0.333) 
C1 = (0.333 X 0.653 + 2 X 0.333 +2 0.653 + 1) 0.766 19.13 x 0.766 | 6 
_ _. (0.653 — 1) + (1 — 0.333) ~ F at 
0:0325. X = 0.043 + 0.33 & 0.0656 | + 0.01 XK 0.333 & 0.653 + 0.257 X 
» 


Determination of Qz: 


2.4nb 


3b 
4 
dy cos ¢ cos ¢ 


G=- + adr > 


when hk = 0, Y, = 0.6: when # = 0.834, Y, = 0.048. 
24X38 xX 100 § (’ x 100 

. = } + 3.14 X 9. 

= + 3.14 X 9.1 X 0.766 | 0.766 


| 425 xX 


= +41.7 (1,682 X* 1,495) — (4,000,000 + 927,000) 


= +49.6 < 108 
(nG)/Ve = (3:8 X 
from Fig. 3, f = 21,400 Ibs. per in.? 


Wy = 


II 


17,500 Ib. 


12 X 100 


0.333 + 0.095 X& 0.653 + 0.105} = $5440. 


+ i) Vs 


y 12/ 
pa — Vin - SP (2 


cos ¢ 


See Fig. 5. 


13) 125 + 12,300 (1 — 0.834)] — 


12 X 100 
0.766 


+ 0.831) o.oss | 


= X 0.60 + 12,300 ( 
0.766 


140,000) “600 = SS7 


(GQ: — Q»)/(f + Q1) = (140,000 x 3,440 — 49.6 X 10°) /(21,400 + 3,440) 
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The Effect of Leading-Edge Bluntness on a 


Laminar Supersonic Boundary Layer 


* 


W. S. BRADFIELD,? D. G. DeECOURSIN,? ano C. B. BLUMER? 


University of Minnesota 


SUMMARY 


An experimental investigation was conducted to ascertain the 
efflect of leading-edge bluntness on the laminar boundary layer 
The results show that the effect of leading-edge 
The ratio of the lead- 
of boundary layer dis- 


of a flat plate 
bluntness is an important consideration 
ing-edge thickness to theoretical values** 
placement thickness at the measuring station will serve as a 
criterion of the importance of the effect It is concluded that the 
presence of a detached shock wave at the plate leading edge is 
responsible for the observed discrepancy between theoretically 
predicted and measured momentum losses 

It is demonstrated that blunting the point of a body of revolu 


tion does not seriously affect the laminar boundary layer 


SYMBOLS 


A, B = constants defined by Eq. (2) 

C, D = constants defined by Eqs. (3) 

Ca. = shock drag coefficient defined by Eq. (7) 

Cy = local skin friction coefficient (t/q .,) 

C; = total skin friction coefficient (f* TwdX/G oo X1) 

h = height of pitot tube tip 

l = distance measured along shock wave 

WV = Mach Number 

p = static pressure 

po’ = impact pressure (measured with total pressure probe 
g = dynamic pressure 

R, = Reynolds Number (p..u.*/p. 

t = flat plate leading-edge thickness or cone tip diameter 
és = temperature 

u = velocity in x direction 


x, Y = coordinates parallel and perpendicular to surface of 


plate, respectively 


8 = Reynolds Number per unit length (R;/x) 

6 = boundary layer thickness (y = 6 at u/m = 0.994) 
5* = displacement thickness [fra — pu/pUq) dy 
6 = momentum thickness 

iv = coefficient of viscosity 

p = density 


Based on a paper presented at the Aerodynamics Session, 
Twenty-First Annual Meeting, IAS, New York, January 26-29, 
1953. Revised and received February 1, 1954. 

* The work reported in this paper was carried out under Con- 
tract Number AF 18 (600) 384, sponsored by Office of Scientific 
Research, Air Research and Development Command 

The writers wish to express their appreciation to Mr. C. J. 
Scott for his cooperation in much of the experimental work and 
his contribution to the data reduction calculations. 

+ Lecturer, Department of Aeronautical Engineering. 

t Associate Scientist, Rosemount Research Center. 

** Throughout the paper reference to theoretical values of 
boundary layer parameters such as displacement thickness, shear 
stress, etc., will mean the idealized theory formulated for a shock- 
less leading edge in the plate case or, correspondingly, an at- 


tached shock wave in the cone flow case. 


T = shear stress 
n = dimensionless boundary layer parameter [yn = 
(y/x)~v ‘R; )] 
Subscripts 
l = conditions at outer edge of boundary layer 
2 = conditions along upper boundary streamline of control 


region 
= conditions in undisturbed freestream ahead of model 


BEL = boundary layer 

m = experimentally measured quantity 

nose = conditions at surface of leading edge of model 

0 = stagnation conditions 

§ = shock wave or conditions behind detached nose shock 
th = theoretical value 

t = conditions at surface of model 


INTRODUCTION 


| Ripepesanegpen INVESTIGATIONS at supersonic speeds 
of total skin friction coefficient by means of the 
momentum loss method have suffered from an evident 
lack of agreement with the laminar boundary layer 
theory as well as a lack of correlation among the results 
of various experimenters. Limitations on the size of 
the models and extent of the laminar boundary layer 
in many facilities have necessitated experimentation 
relatively close to the leading edge of the model. An 
investigation of the effect of flat plate leading-edge 
thickness upon the characteristics of the laminar bound- 
ary layer at M. = 3.05 has revealed a systematic 
dependence of measured momentum defect upon the 
ratio of leading-edge thickness to theoretical displace- 
ment thickness of the boundary layer. 

The present investigation was prompted by the dis- 
covery of the fact that skin friction coefficients obtained 
from momentum loss measurements on a cone exhibited 
good agreement with the theory while poor agreement 
was observed with a flat plate model under the same 
aerodynamic conditions. ! 

Flat plate supersonic boundary layer theory custom- 
arily postulates development of the shear layer from 
a shockless leading edge. Conical supersonic boundary 
layer theory is a direct geometric transformation of 
flat plate theory with the additional consideration of 
an attached shock wave at the cone tip. Actually, for 
blunt leading edges and blunt cone tips a detached 
shock wave must exist in both cases. 

It was reasoned that, due to the three-dimensional 
nature of the phenomenon, the effect of a detached 
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Fic. 1. Flat plate with total pressure probe 
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Leading-edge blunting procedure 


Fic. 2. 


shock wave on the conical boundary layer would prob- 
ably be less than the effect of a corresponding detached 
wave on a plate boundary layer. 

The following experiment was devised to test the 
validity of the foregoing hypothesis. 


EXPERIMENTAL PROCEDURE 


The wind tunnel employed in this investigation was 
an asymmetric channel operating continuously at free- 
stream Mach Number 3.05. The test section was 
13/, in. wide by 2 in. high. The flat plate model was 
made as shown in Fig. 1. The leading edge was hollow 
ground to provide for additional sharpness and to 
avoid the presence of a strong attached shock wave 
on the lower side of the model. The material used was 
tool steel and the model was hardened and ground. 

The investigation included three different leading 
edge thicknesses, namely: 0.3 mil, 5.9 mils, and 11.7 
A sketch of the method of leading-edge blunting 
is shown as Fig. 2. The blunting was accomplished by 
first grinding and then stoning the sharper configura- 
tions to produce the blunter ones. The thickness of 
the leading edge was measured at from 10 to 15 span- 
wise stations in each case. The measurements were 
made by microscope. The thickness figures given 
above are mean values and the maximum deviation 
from these values is +0.1 mil in each case. 


mils. 


AERONAUTICAL 
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This technique of blunting was adopted in order to 
minimize changes in geometric and thermal similarit 
among the configurations of the model. For example 
by preserving plate profile geometry in this manner 
it was expected that alterations in the plate upper sur- 
face temperature distribution due to blunting would be 
minimized. Changes in the geometry of the leading 
edge due to erosion during running were not noticeable 

Thermocouples were installed in a geometrically 
similar soft steel model of 7-mil leading-edge thickness 
at the positions indicated on Fig. 8 and temperatures 
at the upper and lower surfaces were measured. 

Pitot tubes like that shown in Fig. 1 were used for 
the determination of total pressure profiles through the 
boundary layer. Two pitot tube size ranges were used 
in order to permit the assessment of pitot tube size ef- 
fect on the results. The pitot tube tips were rectangu- 
lar in section. The pitot tube tip sizes were 2!,. mils 
by 20 mils, 2*/4 mils by 20 mils, and 8 mils by 20 mils. 
The wall thickness was approximately 1 mil in all cases. 
There was no noticeably different effect of size in the 
data obtained from the 2!/.-mil tube as compared with 
that from the 2*/,-mil tube. 

A photomicrograph of the smallest tube is shown as 
Fig. 5. The response time for the pitot tubes was ap- 
proximately 10 sec.” 

All pitot tube measurements were made at a position 
1*/; in. downstream from the leading edge. Reynolds 
Number was varied by changing stagnation conditions. 
The y-position of the probe in the boundary layer was 
determined with the aid of a 20X microscope during 
running. The possibility of error in probe height de- 
termination due to the presence of density gradients 
and consequent refraction effects was investigated ex- 
perimentally. A very thin flat plate with a scribed 
scale was mounted on the flat plate model at the pitot 
tube position. The orientation of the scale was normal 
to the optical axis of the test section and parallel to the 
Measurements were made with and 


No measurable difference 


flow direction. 
without flow in the channel. 
for the two cases was observed. 

Static pressure measurements were obtained from a 
surface tap in the model surface and by means of a 
static pressure probe in the flow exterior to the bound- 
ary layer. 

Comparable measurements were made on a 15° cone 
in the same channel in the same free-stream Mach 
Number.* 





Fic. 3. 


Photomicrograph of total pressure probe. 
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EFFECT OF LEADING-EDGE 
An investigation of the effect of removing and re- 
placing the model in the wind tunnel on repeatability 
of results was made. The effect on velocity profiles 
was within the experimental scatter of the data. 
Experimental accuracy of the results was as follows: 
pressures +0.015 in. Hg; temperatures +*,, degree; 
viscosity at wall +0.1 percent; pitot tube position +! 
mil average with values ranging from +0.2 mil near 
the surface to +*/, mil at y = 150 mils; Reynolds 
Number +1 per cent average. 
ity profile was found by connecting the first experimen 
tal supersonic point to the origin by a straight line. 


The slope of the veloc- 


Since, in this region, the error in pitot tube position is 
+2 per cent and in u/u, +1 per cent, the corresponding 
uncertainty in the slope is +3 per cent. 

The error in the experimental momentum thickness 
is +3 per cent, consisting of a +2.6 per cent error 
based on measurement of pressure and pitot tube posi- 
tion, and a 0.4 per cent error in the integration of the 
momentum defect through the boundary layer. 

Symbol sizes shown on the plotted results are indica- 
tive of the experimental accuracy in each case. 


DISCUSSION OF RESULTS 


From the total pressure profiles obtained by pitot 
tube measurement, the variation in Mach Number from 
the outer edge of the boundary layer to the wall was de- 
termined. In order to calculate the local velocities in 
the boundary layer, it was then necessary to determine 
the local sound velocities which in turn required a 
knowledge of the static temperature distribution normal 
to the wall. Since no temperatures were measured in 
the interior of the boundary layer, it was necessary to 
make an assumption regarding the temperature distri- 
bution. Isenergic flow was assumed. The effect of 
this assumption on the results was evaluated by com- 
parison with the theory of reference 4 and found to be 
less than | per cent. 


Impact Pressure Profiles and the Effect of Pitot Tube Size 


Fig. 4 is a comparative plot of impact pressure profiles 
with 2' ,-mil and 8-mil tubes in the laminar boundary 
layer on the 0.3-mil leading-edge plate. The effect of 
increasing pitot tube size is to produce an increase in 
impact pressure in the supersonic portion of the bound- 
ary layer to values greater than the free-stream value. 
Hence, a local effective increase in boundary layer 
momentum relative to free-stream values is obtained in 
this region. For thinner boundary layers, the ‘‘super- 
sonic bump” occurred in 2'/.-mil pitot tube measure- 
ments as well as in S-mil A rough 
criterion for avoiding pitot tube size effect in laminar 
plate flow might (according to the present results) be 
tentatively established as h/6* < 0.15. 

The presence of this supersonic bump in the boundary 
layer was first called to the writers’ attention by R. E. 
Blue and G. M. Low whose investigation of pitot tube 


measurements. 
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Fic. 5. Velocity profiles for laminar flat plate flow at 1/, = 
3.05. (Leading-edge thickness = 0.3 mil; pitot tube height = 
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size on laminar boundary layer measurements is pub- 
lished as reference 3. 

The effect on velocity profiles is qualitatively the 
With 


»-mil pitot tube, the supersonic bump is not in 


same as in the case of impact pressure profiles. 
the 2! 
evidence in the laminar boundary layer velocity profiles 
as reference to Figs. 5 and 6 will show. 


Velocity Profiles 


Experimental results obtained from the 0.3-mil lead- 
ing-edge configuration are shown as Fig. 5. By in- 
creasing the ordinates of the theoretical curve‘ by 
12.5 per cent, a curve satisfactorily representing the 
experimental points is obtained. 

It will be noted that the circled points obtained at 
R, = 0.502 X 108 lie below the mean of the three other 
profiles. The fullness of the profile was characteristic 
of all profiles at R, > 0.502 * 10°. As the Reynolds 
Number increased, the profile fullness systematically in- 
creased and the experimental results (see Fig. 7) indi- 
cated transition from laminar flow. The circle points 
were included in Fig. 5 since the value of 26,,/x did not 
However, the effect 
This 


might be interpreted as an indication of an insensitivity 


indicate the onset of transition. 
on local skin friction coefficient is apparent. 


of the total momentum loss calculation (26,,/x) to this 
type and degree of departure of the velocity profile 
from the characteristic laminar shape. 

The effect on velocity profiles of blunting the leading 
edge is shown in comparing Figs. 5 and 6. The bound- 
ary layer thickness is increased as leading-edge blunt- 
ness is increased and the shape of the profile is altered 
near the outer edge of the boundary layer. The curve 
faired through the experimental velocity profiles from 
5.9 mil and 11.7 mil leading-edge configurations was ob- 
tained by increasing the ordinates of the theoretical 
profile* by 36 per cent. 

Increasing the leading-edge thickness from 0.3 mil 
to 5.9 mils results in a decrease in slope of the velocity 
profile (d/dn)(u/u,) near the wall and a departure of 
the experimental points from the adjusted Chapman 
and Rubesin profile at the outer edge of the boundary 
layer. 

A further increase in leading-edge thickness from 5.9 
te ild 
(d/dn) (u/u,) near the wall. 
experimental points far from the wall is marked. 


mils results in no perceptible change of 
However, the effect on 


Close inspection of the 5.9 mil results for 10 < » < 39 
reveals an increasing difference between the experi- 
mental points and the adjusted theory with increasing 
Reynolds Number. This effect is not apparent in the 
case of the 0.3-mil leading-edge configuration. 


Experimental Values of Dimensionless Momentum 
Thickness and Local Shear Stress Coefficient 


Values of 20,,/x corresponding to velocity profiles 
shown in Figs. 5 and 6 are plotted on Fig. 7. It is 
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EFFECT OF LEADING-EDGE 
seen that blunting the leading edge at constant Reyn- 
olds Number causes 26,,/x to increase. 

In evaluating 26,,/x, the departure of measured 
velocity points from the adjusted theoretical velocity 
curve (n < 2) is regarded as pitot tube error. In order 
to account for this effect, velocity values for this por- 
tion of the profile were obtained from a straight line 
joining the first supersonic point of each profile to the 
origin. Values of 26,,/x computed in this manner 
were, on the average, 6 per cent less than values ob- 
tained solely from the experimental points. 

Local shear stress values were computed from the 
relation 

Tw = bw(du/dy) |, 


where du/dy],, was evaluated from the faired profiles 
just mentioned. Values of 4, were calculated from 
surface temperatures measured as described in the sec- 
tion on experimental procedure. 

In the case of the 0.3-mil leading edge, values of C, 
fall below the theory except for the profile represented 
by circled points in Fig. 5 which, as previously re- 
marked, departed from the shape characteristic of all 
the laminar profiles. The onset of transition indicated 
by the local skin friction variation with Reynolds Num- 
ber is consistent with the velocity profile indication 
(Fig. dD). 

The local skin friction coefficient values for the 5.9- 
mil and 11.7-mil leading edges are further below the 
theory and values observed at R, ~ 740,000 for the 
two cases coincide. The latter observation is not sur- 
prising in view of the agreement in the velocity profiles 
shown on Fig. 6 and the thermal similarity of the model 


configurations. 


Heat Transfer Effects 


Due to the shape and material of the plate, heat trans- 
fer from the test surface into the boundary layer is ex- 
pected. The plate was relatively narrow and therefore 
the high recovery temperatures in the turbulent region 
near the tunnel walls together with the high recovery 
temperatures on the lower surface of the plate had the 
effect of raising the surface temperature in the laminar 
The effect of heat transfer on the experimental 
theoretical considerations 


region. 
data was calculated from 
with help from the temperature measurements previ- 
Two cases were investigated: the 
first, described below, assumes isentropic flow at 
the leading edge; the second removes this restriction. 

The leading edge was designed for isentropic decelera- 
tion of the flow along the underside of the plate. As- 
suming that such a flow existed for an insulated two- 
dimensional plate, the temperature distribution along 
the underside would be as shown by the dashed curve 
on Fig. 8. The temperature at the upper surface 
would be the adiabatic recovery temperature corre- 
sponding to the free-stream Mach Number ahead of the 
If the model is of a heat conducting 


ously mentioned. 


leading edge. 
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material, heat must flow from the underside of the plate 
nose to the test surface and thence into the boundary 
layer. In the interest of conservatism (that is, neglect- 
ing conduction effects within the plate) it was assumed 
that the temperature distribution along the lower sur- 
face actually appeared as an equilibrium condition at 
the plate upper surface. The effect of this variation in 
surface temperature on skin friction was then com- 
puted according to the theory of Chapman and Rube- 
sin. The values of Cr,, and C;,, were changed by less 
than 0.1 per cent relative to the insulated plate case. 
The concept of the isentropic leading edge has little 
physical reality since the finest machine work must pro- 
duce an edge which is “‘blunt’’ in the microscopic sense. 
If this ‘‘bluntness’”’ molecular free 
paths of the flowing gas a detached shock wave will 
In this event, temperatures 


is several mean 
form at the leading edge. 
very near the stagnation temperature will occur in 
proximity to the leading edge. 

An estimate of the temperature distribution along 
the test surface in this case is shown as the solid line 
in Fig. 8. Stagnation temperature was assumed at the 
leading edge and the remainder of the temperature dis- 
tribution was estimated from measured surface tem- 
perature as indicated on Fig. 8. 

Thermocouples were installed at the positions indi- 
cated in Fig. 8 on a thermally similar model and 
equilibrium temperatures at the test surface and at the 
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Fic. 9a. Relative difference between observed and_ theo- 
retically predicted momentum thickness as a function of dimen- 
sionless leading-edge thickness from flat-plate measurements 








Fic. 9b. Relative difference between observed and_ theo- 
retically predicted local shear stress coefficient as a function of 
dimensionless leading-edge thickness from flat-plate experiments. 


lower surface were measured. The test points were ob- 
tained from the soft steel model with the 7-mil thick 
leading edge. Reynolds Number per in. was varied 
from 0.187 X 10° to 0.519 XK 108. The effect of Reyn- 
olds Number (hence ¢/6,,*) variation on the tempera- 
ture measurements was less than !/4 per cent. 

The effect on skin friction of the revised temperature 
distribution was computed.* The values of C;,, were 
decreased by 0.7 per cent and values of C,,, were de- 
creased by 0.4 per cent. 

Leading-edge blunting carried out as shown in Fig. 
2 should not alter the temperature distribution along 
the test surface beyond the limits investigated in the 
two examples just described. Therefore, the large 
changes in skin friction coefficients with leading-edge 
blunting cannot be attributed to heat transfer effects. 


Leading-Edge Shock Wave 


The existence of a shock wave at the ‘‘sharp”’ leading 
edge of flat plate models has been established by several 
experimeters.° Fig. 11 is a schematic diagram of an 
assumed leading edge wave configuration. + 
of fluid passing through such a wave would suffer an 


Elements 


initial momentum loss before entering the plate bound- 


{ For convenience in drawing the diagram, the lower half of the 
waves was omitted. 
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ary layer. The total momentum defect of the fluid js 
sensed by the pitot tube at its downstream station 
this includes the initial loss of momentum through the 
shock as well as the momentum subsequently lost io 
viscous shear in the boundary layer. 

The form and location of such a shock wave can be 
expressed to a first approximation solely in terms of the 
freestream Mach Number and the leading edge thick- 
ness.’ A momentum thickness for this wave can be 
defined in a way completely analogous to that for 
boundary layer momentum thickness if the momentum 
loss of all fluid elements which pass through the shock 
prior to entering the boundary layer is considered. 


That is (Fig. 11), 
: ) dl 


which is analogous to the boundary layer momentum 


“5 pu u 
= (1 _ ) dy 
0 Po, U. 


Since 6,, = 0; + @,,, the ratio of 0,/0,, should be a 
criterion of the importance of the effect of shock mo- 


thickness 


mentum losses on total momentum loss measurements. 
Calculation of the momentum thickness for the as- 
sumed wave showed it to be of the same order as the 
theoretical boundary layer momentum thickness for 
conditions of the experiment. 

Since the calculated shock momentum defect is pro- 
portional to leading-edge thickness (for a particular 
M.,)' and theoretical boundary layer momentum thick- 
ness is proportional to displacement thickness, the 
ratio of these two parameters was chosen as a criterion 
for the shock effect in the interest of physical clarity. 

The result of plotting all experimental points from 
the laminar boundary layer against ¢/6,,* is shown in 
Fig. 9. 
tained; for the first, pitot tubes of 2'/. and 2° , mils 


Two complete sets of measurements were ob- 


by 20 mils were employed, and for the second a tube of 
8 mils by 20 mils was used. Over the range of ¢/6,,* 
investigated, the results as plotted in Fig. 9a show an 
effect of pitot tube size on total momentum loss meas- 
urements which is essentially independent of (/6,,* 
variation and which is superimposed on the latter ef- 
fect. The error resulting from the larger ratio of pitot 
tube height to boundary layer thickness has been de- 
fined as negative due to the nature of the distortion of 
the impact pressure and the velocity profiles as pre- 
viously discussed. 

The discussion which follows will be limited to the 
case of the 2!/.- and 2%/,-mil pitot tubes since the effect 
of pitot tube size on the data is thereby minimized. 

The behavior of the curves of Fig. 9 near the origin 
is uncertain. It could be assumed for the purpose of 
discussion that the discrepancy between the idealized 
theory and experiment disappears as ¢/5,,* ~ 0. Ii, 
on the other hand, the curve of Fig. 9a is extrapolated 
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EFFECT OF LEADING-EDGE BL 
linearly to the axis of ordinates, an increment of drag 
exists When ¢/6,,* = 0. This might be interpreted as a 
boundary layer displacement thickness effect at the 
leading edge. 

For ¢ 6y,* > 0.03, the measured momentum loss dis- 
The 
shear stress discrepancy (Fig. 9b) increases rapidly at 
and then less rapidly. Reference 
to Fig. 5 shows that for ¢/6,,* ~ 0.02 the slope of the 
velocity profiles has been decreased relative to the 
theory, but that the shape of the outer portion of the 
profiles conforms to that of the adjusted Chapman and 
Rubesin profile. For 0.32 < ¢/6,,* < 0.56 the slope of 
velocity profiles (Fig. 6) has still further decreased and 
remains constant while the shape of the outer portion 
from the 


crepancy (Fig. 9a) increases linearly with f, 6,,* 


small values of ¢/6,,* 


of the profiles moves systematically away 
theoretical shape with increasing Reynolds Number 
decreasing 6,* for this experiment). At ¢/6,,* = 
1.14, no further effect on the velocity profile slope near 
the wall is seen, but the outer portion of the boundary 
layer has moved quite far from the adjusted pro- 
file. 

Theoretically, the momentum defect (relative to free 
stream momentum) of fluid elements just downstream 
of a detached wave varies comparatively slowly from 
sonic line to sonic line and then decreases quite rapidly, 
approaching zero as the wave approaches a Mach line 
asan asymptote. The modifications of velocity profiles 
just described correspond qualitatively to what might 
be expected if the area of the strong portion of a leading- 
edge shock was systematically increased from zero (by 
progressive leading-edge blunting) while the theoreti- 


UNTNE 


SS ON A BOUNDARY LAYER 379 
cal displacement thickness of the boundary layer was 


held constant. 


Effect of Varying Reynolds Number per Inch Compared 
with Varying Characteristic Length 
If a cross-plot of Fig. 9a is made on a 26,,,x vs. R; 
plot, the results are as shown on Fig. 10. Parameters 
represent constant percentage increments in 
This representation implies that a difference 
x vs. R, curve should be expected 


t/b* 
20 »,/ X- 
in slope of the 26,, 
depending upon whether measurements are made at 
constant with varying Reynolds Number per inch 
or vice versa. Of course, this is important only when 
t 6,,* is sufficiently large that the effect is appreciably 
greater than the experimental scatter of the results. 

if an experimental determination of 


x = 


For example, 
26, x = f(R,, t/6y*) is made at a fixed position on the 
model (x = x, = tf = constant; R,/x, 
variable), ¢ 6,,* will increase as R, increases and an 
increment in slope to the 
will occur as shown by the results of the present experi- 
If, on the other hand, the experiment is set up 
¢ = constant; x; variable), 
and a decrement in 


constant; 


(relative theoretical value) 
ment. 
such that (R, 
t/5»,* will increase as R, 


y = constant; 
decreases 
slope should be expected. Computation from the ex- 
perimental data (5.9 mil leading edge) of the expected 
variation of 26,,,x with R, under these conditions yields 
the curve shown on Fig. 10. Qualitative confirmation 
of this effect is shown by the experimental results of 
Maydew and Pappas.* 

In view of the dependence of 26,,/x on the manner 
in which Reynolds Number is varied, it is of interest to 
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DEFINITION OF CONTROL REGION FOR MOMENTUM LOSS MEASUREMENTS 





Me «3 
te 
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poe STATION, X, 
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— if =X = on 
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fe) 0. 0.2 0.3 0.4 0.5 
Ee | 1 J 
SCALE. INCHES 
Fic. 11 Definition of control region used in momentum loss measurements 
extend this investigation to the case of local skin fric- Cy, = —2[(B/x)(00,,/OB) | 


tion coefficient determination from momentum loss 
data, 


As long as the pitot tube is sufficiently far from the 
leading edge and the plate is aligned with the free- 
stream flow, dp/dx| = 0. Therefore, 

= X) Cy = 2[(d0»\/dx] 
and it should be possible to compute local skin friction 
coefficients equally well from the equivalent expressions 
Cy = 2[(d0,)/dx] = 1/q , \ou(du/dy) |r} 

Figs. 7 and 10 show that where (/6,,* 

(for a particular 7 .), 


Om = f(B, x, t) (1) 


is appreciable 


Referring to Fig. 9a, the linear portion of the curve 
may be described by 
Om = 6y,{At, én" + (B + 1)] (2) 

where A is the slope and B is the intercept of the ex- 
perimental data. 

Oi, = CB '<; . 

bn* = Dpo "x" (3) 
for a particular free-stream Mach Number. Substitut- 
ing and differentiating, 


(OO m, Ox) Jp. ,_= | 2(B + 1 CB 1 - 


B a | es 
= —(1/2)(B + 1)CB~"“x 
* OB |. 4 
That is 
Cy = 2(06m/Ox) lg. 4 = —2}(B/x)(0m/O8)} |x, (4) 

Eq. (4) holds as long as Eq. (2) is linear in (¢/6,,*). If 
it is not linear (which may be the case for 0 < t/6y,* < 
0.03), Eq. (4) no longer holds and an expression which 
is explicit in ¢/6,,* will result. 

Since the form of this curve is not completely defined 


by the present experimental results, a linear relation 
was fitted for ¢/é»,* > 0.03 and values of 


were obtained from the experimental results. De- 
pending on the fairing of the experimental points plot- 
ted as 6,, vs. 8, local skin friction coefficients were ob- 
tained ranging from 40 per cent below the theory to 
10 per cent above the theory. The sensitivity of the 
calculation to fairing is such that unless the accuracy of 
the experimental values of @,,/x and R, is extreme or a 
large number of experimental values is obtained, the 
technique is inaccurate for determining local skin fric- 
tion coefficients. 


Evaluation of Shock Drag Coefficient 


If a leading-edge shock wave is postulated as the ex- 
planation for the blunting effect, the region of interest 
in the flow would appear as shown in the scale drawing 
of Fig. 11. In flat plate boundary layer experiments, 
it is commonly assumed that the total defect of momen- 
tum through the boundary layer is due solely to the 
action of the shear forces at the surface. In the present 
case, the difference between the total momentum loss 
determined from experimental points and theoretically 
predicted values of total momentum loss due to skin 
friction would include the shock drag. That is, 
2(Om/X — O1,/X) 

Crm 
were drag coefficient + ~~ C 
_ F th 


| total skin friction coefficient 5) 
(oO 


on 
Referring to Fig. 11, the total change through the 


control region of x-momentum as calculated from the 
pitot tube measurements is related to the forces by 


20» 2 5+ 1/2 pu u 
= 1 — dy 
v1 xX t/2 pu. u 
l “1 1 
=— Cy(x) dx + -~ x 
x1 Jo qo¥1 


t/2 : l 5 + 1/2 
f (Prose — p a dy = [ (po — Po) dy 
0 G.%1 Jt—* 


(6) 
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EFFECT OF LEADING-EDGE BLI 
where 
p pressure along the upper boundary streamline 
of the control region 
q l 2p i 


The last term on the right-hand side of Eq. (6) repre- 
sents the pressure force due to possible pressure gradi 
ents in the flow exterior to the boundary layer, down- 
stream from the shock wave. The conditions under 
which this term may become important will be dis- 
cussed in what follows. The control region .properly 
defined is bounded on the lower side by the stagnation 
streamline and the test surface and on the upper side 
by the first streamline conforming to the condition that 
V M., pi ~ p 


The shock drag coefficient is given by 
5 5 » 


» 2ee. 


where 


is the true total skin friction coefficient. 

If the physical picture just described is correct, it 
should be possible to compute the shock drag coefficient 
from Eq. (7). The last term of Eq. (7) will not enter 
the computation if the traverse at x = x 1s extended 
sufficiently far in the y-direction to ensure that the 
upper boundary streamline of the control region passes 
through the Mach wave asymptote of the leading-edge 
po =~ p, and the 
However, the 


shock. Under this circumstance, 
contribution of the term is negligible. 
asymptotic character of the impact pressure profile at 
x = x; makes it difficult to define the outer boundary 
exactly and small discrepancies between local impact 
pressure and free-stream impact pressure may result 
in values of pop. appreciably greater than unity. 

For example, in the case of the 11.7-mil leading edge 
the traverse at x = x, was extended out to 145.5 mils. 
This compares with 6 = 22 mils at the same R, with a 
0).3-mil leading edge. In constructing the upper bound- 
ing streamline for this control region, it was found that 
it passed through a portion of the leading-edge wave 
which was of finite strength. The ratio of impact 
pressures through this section of the wave was Po,/ Po 
= 0.998. However, the computed ratio of static pres- 
sures was p,/ p,, = 1.5. 

The discussion of the last two paragraphs is based on 
the assumption of a hyperbolic shock wave configura- 
Therefore, the quantitative value of the discus- 
However, the 


tion. 
sion must be tempered by this fact. 
comparison of relative changes in impact and static 
pressures through weak oblique waves indicates the 
qualitative importance of this consideration. 
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Summary Discussion 


7) from the experimental re 
the first two 


An evaluation of Eq. 
sults is shown in Fig. 12. The sum of 


terms is shown as the solid curve. The total skin fric 


' a wg 
Cr = / CrAx) dx 
v} 0 


. 


tion cé efficient 


was computed from the experimental values of C, 
plotted on Fig. 9b. A hyperbolic variation of C, with 
t 6,,* was assumed in fairing a curve through the ex 
perimental points. 

The pressure force on the upper boundary of the con 
trol region (the third term of Eq. (7)) was estimated. 
The initial value of p2 was obtained as described above 
The integral was evaluated assuming a linear decrease 
from the initial value of po = p, to the measured values 
of fi = p, atx = M1. The bars (I) indicate the ap 
proximate nature of the calculation, and the length of 
the bar represents the maximum expected contribution 
of this term in each case. 

Shown for comparison with the experimental results 
is the shock drag coefficient calculated from the theory 
of Moeckel’ which has been proved satisfactory for 
blunt two-dimensional and axially symmetric bodies 
in an inviscid flow. The agreement with the present 
results is fairly good considering the uncertainties in 
evaluating Eq. (7) from the experimental results. 

For example, the definition of the control region for 
the experiment may have a large influence on the evalu 


ation of Eq. (7). In addition to the previously dis 
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cussed possibility of an undetected force on the control 
region upper boundary, there seems little doubt that 
the position of the leading edge stagnation line is an 
important consideration. Migration of the stagnation 
line to the upper or lower surface of the model would 
probably affect the contribution of the shock wave to 
the total momentum loss along the test surface. 

Furthermore, the calculated value of shock drag coef- 
ficient shown for comparison with the experimental 
results is proportional to f/x,. The experimental ac- 
curacy of leading-edge thickness was +0.1 mil and a 
corresponding uncertainty exists in Cua,,,. Finally, 
Moeckel’s analysis was not intended for cases such as 
the present one where viscosity may play an important 
role in the leading-edge vicinity. A more suitable 
theory does not exist, at present, to the writers’ knowl- 
edge. 

Correction of experimental momentum loss data to 
obtain values of Cy and C, for comparison with the 
idealized theory would imply a reversal of the proce- 
dure required to obtain calculated values of C;,. Hence, 
it appears that calculation of Cry and C, data from 
momentum loss experiments by a correction procedure 
would be of doubtful value. 

Differentiation of measured total values of momen- 
tum loss to obtain local shear stress coefficients was too 
inaccurate to produce any conclusive information from 
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this experiment. However evaluation of local values 
of viscous shear at the plate surface was obtained by a 
consistent estimation of the velocity profile slope at the 
wall together with values for viscosity from measured 
temperatures as previously described. As shown by 
the results of Figs. 9b and 13b local shear stress coef- 
ficients so determined are not sensitive to changes in 
t 6,,* for values of ¢ 6,,* greater than about 0.3. 

For a pointed body of revolution, the comparative 
effect of a detached tip shock on the laminar boundary 
layer is much smaller than is the case with flat plates. 
This is explained by the difference in geometry of the 
flows and of the shock configurations. The corre- 
sponding calculated ratio of Cy,/Cr,, for a 15° conical 
body at J/, = 3.05 is shown for comparison with the 
plate calculations in Fig. 12. A summary of experi- 
mental results of the cone blunting investigation is 


shown in Fig. 13. 
CONCLUSIONS 


The results of the present experiment indicate that 
the effect of leading-edge bluntness is an important 
consideration in the formation and development of the 
laminar boundary layer of a flat plate at supersonic 
speeds. 

The effect of the leading-edge shock wave is to cause 
a thickening of the boundary layer and to reduce the 
slope of the velocity profiles adjacent to the wetted 
surface relative to what would be expected in the un- 
disturbed or shockless case. The former effect ex- 
hibits a leading-edge 
thickness for fixed aerodynamic conditions. The lat- 
ter effect, however, exhibits no discernible systematic 


systematic dependence upon 


dependence upon the variations of leading edge thick- 
ness. 

The ratio of leading-edge thickness to theoretical 
values of boundary layer displacement thickness at the 
measuring station may be taken as a criterion of the 
magnitude of this effect. 

The blunting effect is shown to be separate from the 
effects of heat transfer and pitot tube size. 

It is noted that the effect of varying leading-edge 
bluntness on total momentum loss in the boundary 
layer is generally stronger than the effect on local 
values of viscous shear at the plate surface (7, = 
u,du,dy),,) estimated as described above. 

It is concluded that the presence of a detached shock 
wave at the plate leading edge is responsible for the ob- 
served discrepancy between the momentum loss pre- 
dicted by the idealized theory and measured momen- 
tum losses. 

Due to the experimental difficulty in accurately de- 
fining the pressure drag on the control region, it 1s 
doubtful that correction of observed data for the shock 
wave effect would be successful. In order to obtain 
correlation between the plate boundary layer data and 
the idealized theory, it is necessary to minimize ¢ 6,,* 


(Continued on page 398) 
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Application of Volterra Linear Integral 
Equations to the Numerical Solution of 
Vibration Problems—IT 


J. L. BOGDANOFF,? J. E. GOLDBERG,? ann HSU LO*™ 


Purdue Unwersity 


ABSTRACT 


A new and simple numerical method is presented for the solu- 


tion of one-dimensional characteristic value problems. Three 


examples dealing with beam vibration are given 


(1) INTRODUCTION 


on NUMERICAL RESULTS GIVEN in the present paper 
are a further application of a method for the ap- 
proximate solution of characteristic (or boundary) 
value problems developed by the authors.’ Very 
briefly, the method consists in replacing the differential 
equation originally present in the problem by a linear 
Volterra integral equation, and solving the integral 
equation numerically subject to the boundary condi- 
tions. In the previous paper, only problems involving 
linear differential equations with constant coefficients 
were considered. The examples considered herein have 
variable coefficients in the differential equations and 
pertain to the torsional and lateral vibrations of thin 
beams of variable section. The paper begins with a 
brief discussion of the basic relations used in the ap- 
proximate method. 


(II) Baste RELATIONS 


Consider, for example, a characteristic value prob- 
lem, defined by the linear differential equation, 
(d4y/dx*) + a3(x)(d'y/dx*) + ....4 


Ady(x)y = O, oss 5 


and the homogeneous boundary conditions 


y(O) = QO, (=) = 0| 


Gx «a0 
d*y d*y 
3 = {), - a= () 
Gi" 2 a% Gn! os 
where a;(x 0, 1,...,° are known functions of x, 
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and \ is a parameter. (With proper interpretation of 

the symbols, Eqs. (1) and (2) could be used to determine 

the natural frequencies and principal modes of lateral 

vibration of a thin nonuniform cantilever beam. ) 
Introduce the substitutions 


dty/‘dx* = Y(x) 
x 
a’y dx* = f V(é) dé — C, 
0 
d*y dx? = / (x — &) V(é) dé — Cw ( 
70 
"= (x — §)* x (3) 
dy dx -f aa V(é) d— — C5 - Con 
0 - 
“x(x — &£) 
vv) = : V(t) dé — 
./70 6 
x° _¥ 
C; " cad C2 e 
6 2 


These rela- 
The 


where C, and C, are arbitrary constants. 


tions clearly satisfy the first two of Eqs. (2). 


Volterra linear integral equation of the second 
kind 
V(x) = Cif(x,A) + Cog(x,A) — 

x 

K(x, &; A) ¥(&) dé (4) 

J 0 
where 
f(x,rA) = as3(x) + xa2(x) + (x?/2)ai(x) + 
A(x? /6)ao(x) 

g(x,A) = do(x) + xay(x) + A(x*/2) ao(x) > {> 


a3(x) + (x — E)ao(x) + 
2]a,;(x) + A[(x — &)8 


K(x,&; A) = 


[(x — £)? 


6 lao(x) | 


results from substituting Eqs. (3) into (1). The third 


and fourth equations of Eqs. (2) now become 


1 
f (1 — &) V(t) di —-C, -G= 0 
ys (6) 


. -_ 


| V(é) dé — CG, = O 


The characteristic value problem originally defined by 
Eqs. (1) and (2) is now formulated in terms of the Vol- 


terra integral Eq. (4) and Eqs. (6). 
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The trial-and-error method used by the authors to 
obtain approximations for the first m-characteristic 
values and functions is the following: 

(1) Divide the interval (0,1) into n-equal parts each 
of length h. 

(2) Select a value of \ for trial; denote it by A*. 


(3) Compute the f(x, A*) and g(x, A*) at x = 0, 
h,....,th,...,nh = 1, and denote them by 
er | F,, 
nis <3 net Ce 


respectively. 

(4) Let €. = ft and C,:=.0: 
follows from Eqs. (4). 

(5) The value of V(x) at x = h, V4, is determined by 


Eq. (4): 
*h 
y, = Ff, - / K(h,é; X*) V(E) dé 
0 


An approximate value for }; can be found from this 


Then, Y (0) = F,; this 


equation by assuming ](£) varies linearly between 0 
and h, i.e., 


V(é) ~ Yo + [(%1 — Vo)/A]é, O<ét<h 


evaluating the integral with this assumption, and then 
solving the resulting algebraic equation 


y; to F, — A 10 Yo = Ai VY; (7) 
for an unknown J;: 
Y, = [(F; = AwYo) (1 + Ay) ] (8) 


The known quantities Ay, Ay are defined in Appendix 
I. 

(6) To find an approximate value for V(x) at x = 
2h, Eq. (4) is again used: 


2h 
Y, = Fy - | K(2h,€; \*) V() dé 
0 


The integral may be approximately evaluated by as- 
suming 


V(é) ~ Yo+ [(%1 — YVo)/hlé, O<é<h, 
V(é) ~ Vi + [(¥2 — Vi)/h] (& — h), O< E < Qh. 


¢ 
Then, the above equation becomes 


Y, —_ Fy —_ AnYo = An Vi = AwY> 


l 
Y2 = (¥: as ® Au.) [( + A») (9) 
1=0 


Continuing in this manner, always assuming that Y(&) 
may be approximated by means of a piecewise linear 
function in evaluating integrals, one can find the ap- 


proximate value of V(x) at x = 3h, 4h,..., nh = 1. 
The general formula is 


r—1l 
Y,=~ (F. — DAs, v,) [a + A,,) (10) 


i=0 


and 
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Denote this sequence by 


SPS Oran 


(see Appendix I). 


(7) Now let C; = 0, C2 = 1. Then Eq. (4) shows that 


Y (0) = G (0) 
(S) Find the approximate value of V(x) at x h, 


2h,..., nh = 1, by the procedure outlined in steps 5 
and 6, using the general formula 


, ] 
¥.~(G - Daur) (1 + A,,) 11) 
1 0 
Denote this sequence by 
Dis Bee sone 
(9) Then, the general Y(x) at the division points are 
Y.~ GY,’ + CY’, i=0,1,...,n” (12 


Evaluate approximately the integrals appearing in 
Fqs. (6) using these Y; and assuming Y(£) varies lin- 
early between adjacent division points. The substi- 


tution of these values in Eqs. (6) gives 


PiC; + PoC, = Ol 
Q:C; + QoC, = OF 


(13) 


where P;, Ps, Q:, Qe are defined in Appendix IT. 

\* is a characteristic value (approximately) if the de- 
terminant of the coefficients of C; and C. vanishes. As 
a rule, the initial value selected for \* will not cause 
this determinant to vanish. However, after a few trials, 
a value can be found which is satisfactory. 

(10) Let A, denote the approximate value of the kth 
characteristic value and let the corresponding values 
of VY; be denoted by 


Vin = GV ux’ + CVx”, p= 6, 1,...,0 1 
Since the determinant of the coefficients in 

Pi.C; + PoC. = 0 (15 

”) 


Oud: + OxC: = Of 
vanishes 
C,/C2 = —(Px/Px) = —(Qx/Qx) 
Thus, 
Vin = Cit Vue’ — [(Pic/Pax) Vie" VI, ¢= 0, 1,...8 


(11) The approximate values of the corresponding 


characteristic function at x = 0,h,..., nh are found 


from the last of Eq. (4): 


F th (th — £)% Y(é) th)3 
ys _ a=2 ( az - (th 4 
Cy 0 6 C; 6 


P,. (ih)? 


P tae . 2 


Again, the integral is approximated by assuming that 
Y(é) varies linearly between adjacent division points 
and at the division points equals 
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f(x,A) = Ax + ay(x), 





APPLICATION OF VOLTERRA 


Pen i) ge 


This, in very brief form, is an outline of the method 
used by the authors to obtain the characteristic values 
A more de- 
Clearly, it 
was not necessary to consider a fourth order linear dif- 


and functions of the systems (4) and (6). 
tailed treatment is given in reference 1. 


ferential Eq. (1) and the particular homogeneous bound- 


ary conditions, Eqs. (2). 


(III) EXAMPLES 
(A) Torsional Vibration of a Fixed-Free Beam of 
Variable Section 
The physical properties of the rod considered in this 
example are shown in Fig. la. The characteristic value 
problem for this example is defined by 


(d?0/dx?) + [a\(x) (dO/dx)|] + 40 = 0, <<. f, 
(¢ ? 
6(0) = 0, = 0 
dn 

where 

Z = coordinate along axis of rod 

L = length of rod 

x = 3/L 

6(x,t) = O(x) cos pl 

a(x) = —5.544x 

a(x) = 1 

p = mass density 

G = shear modulus 

d = (plL*p*)/G 

Vix) = d*0/dx? 

C(x) = torsional stiffness 

I(x) = mass moment of inertia 


The equations analogous to Eq. (4) and Eqs. (6) 


are 


Y(x) = Cif(mA) — i K (x,&; A) V(&) dé 
0 


wa 
0 


respectively, where 


K(x,€; A) = 
A(x — £) + a(x) 


Approximate values for the first four angular natural 
frequencies and corresponding mode shape values at 
the division points, with » = 10, are given in Figs. 1b, 
le, Id, and le. The fifth frequency and mode shape 
with x = 10, = 20, and = 40 are given in Fig. 2. 

The fifth frequency is accurate to within about seven 
per cent with x = 10, the lower order frequencies being 
more accurately determined. Fig. 2 also reveals that 
with » = 10, the node points are quite accurately deter- 
mined. 


L 
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Fic. 1. Torsional frequencies and modes of tapered fixed-free 
rod; 10 stations 
(B) Torsional Vibration of a Fixed-Fixed Beam Carrying 


a Rigid Disc 
To demonstrate that the method works equally well 
with discontinuities in the a;(x), consider the system 
illustrated in Fig. 3a. The equations defining the char- 
acteristic value problem are: 





(d?0 /dx*) + [a;(x) (d6/dx)] + AO(x) = 0, } 
0s. 25 1] (a) 
(d*@/dn*) + AP(n) = 0, O<7< 1 \ 
and 
6o(0) = O, 0(1) = (0) 
1m I, Lp? 10 
( ) = ~~ aa) = ( ), oa) = 0 (b) 
dn 0 GJ» Ix 1 
where 
\ | 5 Move 


10 STATIONS —~ 


of 4 ; i x 
ar) a 
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Influence of station number on 5th torsional frequency 
and mode of tapered fixed-free rod. 
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d = (pL*p?)/G 
O(x) = angular amplitude of cross section of rod to 
left of rigid disc 
a(x) = —2.7726x 
#(y) = angular amplitude of cross section of rod to 


right of rigid disc 


and assuming that the shafts have circular cross sec- 
tions, are of the same material, have equal lengths, and 
the diameters are the same at the rigid disc. The in 
tegral equations for (a) are 









































V(x) = CitAx — 2.7726x} — 
wy 
J tA(x — &) — 2.7726x} Y(é&) dé 
0 (c) 
a) 
Z(n) = DirAn + AD, — / Ain — £)Z(¢) de 
0 
where }(x) = d°O/dx’, Z(n) = d*/dn? and the first 
of (b) has been used. The second and third of (b) 
show that 
D (=) » © a1) 
>= sail 7 
dx/ pLJy (d) 
Finally, the last of () is equivalent to 
Zl) = 0 
because of the second of (a). 

In the computations, /)/ pLJ) was taken equal to 
unity, and # equal to 1/5. Results for the first two fre- 
quencies and mode shapes are shown in Fig. 3b. 

(C) Lateral Vibration of a Nonuniform Cantilever Beam 

The equations 
a de ae ae ae 

oy, ili =? ne aay = 0, rs a) 
dx? | dx? § 
and 

oo : F 
r tk, ' 
a) Jeurqt= Ge" 4 GI, = const b 
a 4eT, <7 wens xt pJe= const f 
+> x e n 
t 
role 
° 
3 . ' 
° | ‘ 
b) ee 
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° UE L 
«x \ 
° % ! 
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Fic. 3. Torsion frequencies and modes of tapered fixed-fixed rod 


with rigid disc at center; 10 stations. 


teCeatL 
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W(0) = 0, 


7) 
= (0), 
fe ae 


define the characteristic value problem for the example 
under consideration. The definitions of the new sym 


bols are as follows: 


w(x,t) = lateral displacement of centerline 
w(x,t) = W(x) cos pt 

jm = pAoL'p?/(EI), 

Ay = area of cross section at x = 0 


(EI)o = bending stiffness at x = 0 


With substitutions of the type given in Eqs. (3), (a) 
and (b) become 


Q(x) = Ci {2cex(2ex? — 1) — 4ex — p(x?/6)} + 

Cz {2c(2cx? — 1) — p(x?/2)} — 
x j 
f 42c(2cx? — 1) (vw — &) — 4exv - 
0 
(x — &)%) 
Me > Q(—) dé (ce) 
6 f 
and 


"" 
J (x — )0(¢) dé - CG. -G = 0] 

0 

- (d 
| Q(t) dé — C, = 0 | 

0 


respectively, where Q(x) = d‘*W/dx! 

Results obtained for first four angular natural fre- 
quencies in lateral vibration and mode shapes, with 
10 and c = The first 
and third frequencies and modes shown in Fig. 5 refer 


n= 1.25276, are shown in Fig. 4. 


sO 


to the case in which c = 2.772. 


(IV) DiscussIon 


There are several points of interest concerning the 
method. First, the Volterra integral equation is re- 


To 


placed by an approximately equivalent system of 
linear simultaneous nonhomogeneous algebraic equa- 
Since the 
matrix of the coefficients of these unknowns is triangu- 


tions in the m + 1 unknowns Vo, i,..., Y,. 


lar in form (see Appendix I), the equations are easily 
solved progressively. Second, it follows from the forms 
of the system of linear equations and boundary condi- 
tions, and the fact that only n-subdivisions are used, 
that estimates can be obtained for at most only the 
first Third, 
the order in which these estimates are obtained is im- 


n-characteristic values and functions. 


material. Thus, if it is desired to compute the Ath 
(<n) characteristic value, it is not necessary to de- 
termine either the first k-1 characteristic values or func- 
tions. It follows that the labor required to obtain k 
characteristic values and functions is roughly propor- 


tional to k. Fourth, the boundary conditions and not 
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APPLICATION OF VOLTERRA 

the number of subdivisions of the interval (0,1) control 
the order of determinants which must be used in the 
problem. The only determinant which it is necessary 
to expand in section II occurs in connection with Eq. 
(13), and the order of this determinant depends only on 
the number of boundary conditions. This is important 
when using a desk calculator and an unskilled com- 
puter, and is in marked contrast to certain other nu- 
merical methods which start from a Fredholm integral 
equation formulation of the characteristic value prob- 
lem and replace the integral equation by an approxi- 
mate system of linear simultaneous homogeneous alge- 
braic equations. Fifth, there are two distinct approxi- 
mations involved: both the integral equation and the 
boundary conditions at x = | are satisfied only approxi- 
mately. In spite of these two approximations, experi- 
ence seems to indicate that the estimated characteristic 
values are always greater than the exact values. The 
reason for this has yet to be assigned. Sixth, the 
method is simple to use with a desk calculator. For 
once tables have been constructed for the basic elements 
appearing in the A,; (and this has been done), it is easy 
to set up a tabular procedure for carrying out the indi- 
The method is somewhat reminiscent 
Seventh, 


vidual trials. 
of, though not related to, Holzer’s method. 
the quantity A(x,£;) occurring in Eq. (4) is obtained 
in a routine and simple manner. If Eqs. (1) and 
(2) had been replaced by a Fredholm linear integral 
equation, it would have been necessary to determine 
the Green's function for the system; this is occasionally 
a matter of some difficulty. Eighth, the method is ap- 
plicable if the a;(x) are given in graphic or tabular form 
and if they are piecewise continuous. Finally, the 
method can also be used to determine the particular 
solution arising from a constant nonhomogeneous term 
in the boundary conditions with X specified. This 
follows directly from the forms of (4) and (6). It also 
implies utility of the method in forced vibrations prob- 
lems where a known excitation (simple harmonic) 
acts at a boundary point. 

After using this method for over a year, the authors 
found that in 1939 A. Huber? used the same type of ap- 
proximation for }(£) in obtaining an approximate nu- 
merical solution of a Volterra integral equation of the 
first kind which arose in connection with a heat conduc- 
tion problem. His results were excellent. Huber also 
made some comments on the nature of the error in the 
method which do not appear to apply under general 
conditions. No comment was made, however, on the 
utility of the method in characteristic value problems. 

L. Collatz* has suggested that an estimate for the 
first characteristic value can be obtained easily by means 
of a successive approximation method that starts from a 
Volterra integral equation formulation of the character- 
istic-value problem. Although this suggestion does 
not appear to be well known, it provides a useful alter- 
native to the well-known Rayleigh method. Clearly 
Eq. (4) could be solved, subject to Eq. (6), by successive 
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approximation. The authors did not pursue this 
method, however. 

It should also be pointed out that it is not necessary 
to use trial-and-error in finding the A,.__ Instead, A could 
be carried along as a parameter in evaluating the Y 
Then, the vanishing of the determinant of the coeffi- 


cients of the C; and C, in Step 9 would yield an algebraic 
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equation of mth degree in \; and the roots of this equa- 
, A, The 


corresponding characteristic functions could be found 


tion then provide approximations for \,, . . . 
as outlined in Steps 10 and 11. The authors tried this 
approach initially, but gave it up because of the heavy 
algebraic labor involved. 


V. COMMENTS ON ACCURACY 


Experience indicates (see example A and reference 
1) that engineering accuracy can be obtained with 
roughly twice as many stations as the order number 
of the highest characteristic value wanted. Suppose, 
in order to be definite, that it is necessary to determine 
the fourth natural frequency of a simply supported beam 
in lateral vibration. To secure engineering accuracy, 
eight to ten division points should be used. If extreme 
accuracy is required, the authors suggest that the fre- 
quency be computed twice, once with eight divisions 
and once with twelve divisions, and then use the method 
of the deferred approach to the limit to improve the 
results. At least four divisions should be used in esti- 
mating the first characteristic value. Examination of 
the examples reveals that the characteristic functions 
are unusually well determined. In particular, the 
nodal points are accurately determined even with a 
small number of divisions. This is contrary to results 
obtained by other methods. As is to be expected, the 
error increases as the order of the characteristic number 
increases. 

The use of a piecewise linear approximation for }’ (€) 
in evaluating integrals is clearly not the most accurate 
approximation possible, although it is quite simple to 
use; and numerous schemes for improvement suggest 
themselves. For example, it is possible to use a para- 
bola for Y(é) in the interval 0 < — < h to find Yj, ad- 
justing the slope at the origin so that it equals dV ‘dx at 
x = Owith C, = 1, C, = 0, and having the parabola pass 
through VY) = Fy. The function F(~, \*) could also 
be used for }¥(£) in the first interval. 
a cubic parabola through }y, 1, and with the proper 


In either case, 


slope at x = 0, would permit the determination of Yo, 
etc. Once tables for the basic elements in the A,; had 
been constructed, the procedure would not be any more 
complicated than that suggested above. Alternatively, 
one might use some of the methods suggested by G. 
Prasad,‘ and L. Fox and E. T. Goodwin.’ Nothing 


has been done in this direction, as yet. 


APPENDIX I 


If in step 5 of section II, the linear approximation is 
substituted in the integral equation, the equation 


h , 
Yy~F, - f K(h,£; X*) (% a Lf j Vo e) ae (a) 
0 H 


results. Hence 


Y, a F, —= AwYo as An’V; (7’) 


TLCAL 


SCIENCES JUNE, 1954 


Ay = J Khe: \*) (: ~ 4 de] 
0 h 

A, = | K(h,&; \*) . dé | 
0 h 


Since these A’s are known quantities and Vy is also 


where 


known, Eq. (7’) can be solved for Y; (see Eq. 8). 
Suppose that Yo, ,:.., VY,; are known and it js 
desired to find Y,. From Eq. (4) 


yy, = F, — /  K(rh,t: A*) V(&) dé (e 


The replacement in this equation of Y(é) by a piece- 
wise linear approximation gives 


where 


Aw» 


‘ | 
A, = ft. K(rh,é; *) ( —it+ i) dé + 
t—Il)h h 
_ 


K(rh,&; X*) (: 


ctr 


2 . c (e) 
K (rh,&; X*) (i +] — *) dé, 


Str 


*rh 
Ar = i K(rh,é; \*) ( —r+ 1) dé 
) lh h 


Eqs. (9) and (10) follow directly from (d). Once 


K (x,£; \*) is known and X* and h selected, the A,; can 


be computed. The general form of (d) or Eq. (9) is 


Yo = Fo 

Aw¥Yo t+ (1 + An) VY = F, 

bhikt hee + ...%64+40%=% 

Ae¥o + Aa¥i +.....-+0+4+40)¥, =F, 
Since the matrix of the coefficients is triangular in form, 


these equations can be solved directly without having to 
resort to the use of determinants. 


APPENDIX II 


» 


The substitutions used to arrive at Eq. (13) show 
that 


h ra 
P, = vf (1 — £) (1 — s) dé + rs x 
0 h 
f h é “2h :) ) 
1 — §)- dé (1 —€&(2- ltée + 
if ay + | ( nl) °*§ 


! é 
a—9( —n+1)ae-1 
n—I)h h 


(Continued on page 403) 
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Extreme Values in Aeronautics 


EMIL J. GUMBEL* ano PHILLIP G. CARLSON? 
Columbia University and Curtiss-Wright Corporation, Respectively 


SUMMARY 


It is known that large and small gust accelerations encountered 
in flights at essentially constant airspeed are distributed according 


to exponential laws. Hence, extreme accelerations, gust veloci- 


ties and maximum airspeeds observed under certain conditions 
may be analyzed by the asymptotic theory of extreme values 
To this end, a short résumé of this theory, graphical procedures 
for its use and for the control curves are given in the first para- 


graph. In paragraphs 2 and 3 it is shown that extreme winds, 


maximum airspeeds, and extreme gust accelerations follow this 


distribution. Finally, the theory is applied to largest, smallest, 


and extreme accelerations. 


INTRODUCTION 


a THE ENGINEERING POINT OF VIEW the theory of 
extreme values is most useful for design. The civil 
engineer planning a river dam or bridge must estimate 
the largest floods which may occur during the expected 
life of the structure. The aeronautical engineer is 
interested in, among other things, the probability of an 
extreme load which might cause structural failure. He 
is concerned with the probability of gusts exceeding a 
prescribed magnitude over a given period of time. His 
are the problems of finding the effects of the frequency 
of large loads on the operational life of an aircraft, or 
the probability that a large gust will cause unstable 
operation of the power plant of a supersonic aircraft. 

On the other hand the engineer finds that the smallest 
value leads to fruitful applications in the statistical 
analysis of droughts and the breaking strength of 
materials. 

These applications of the theory of extreme values to 
aeronautical problems are by no means exhaustive. 
Y. C. Fung* has shown that the distribution of extreme 
values can be used to construct an “‘envelope’’ repre- 
senting the most probable largest stress and other 
response quantities of interest at any point in the 
structure which, used as a basis for specifying the design 
If the 
probability of exceeding a given maximum stress is 
known, the life expectancy of an aircraft can be deter- 
mined as a function of the specific load and vice versa. 


stresses, leads to a uniform factor of safety. 


(1) The Theory of Extreme Values 

Among » observations of a statistical variate, there is 
one, x;, Which is the smallest value, and another, x,, 
which is the largest. In a second series consisting again 
of » observations made under the same conditions, the 
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smallest one will be, say, x,’ and the largest x,’. If NV 
series are observed, .V smallest and .V largest values are 
obtained. These are called the extremes. The problem 
is to find the distributions of the extremes, especially 
if is large. 

Let F(x) be the initial probability function. 
the probability that » independent observations are 
below x is the product F"(x), and the probability that 


Then 


all observations exceed x is the product || — F(x) 

The first expression can be interpreted as the probability 
®,(x) that x is the largest among m observations. The 
that x 
The derivatives 


second expression is the probability 1 — \®,(x 
is the smallest among m observations 


¢rn(x) = nF"—'Mx)f(x), re n(x) = n{l — F(x) ]"—' f(x 
are called the distributions of the largest value x X 
and the smallest value x = x,;. The probability func- 


tions ®,(x) and ,®,(x), and distributions ¢,(x) and 
1¢,(x) depend upon the sample size » and the initial 
distribution f(x) = F’(x) from which the samples are 
taken. 

For large values of the sample size m there exist three 
different asymptotic expressions.‘ Considered here is 
the case where the initial variate is unlimited in the 
direction which is of interest. For large positive values 
of x it is assumed that the initial probability function 
F(x) approaches unity in the same way as e~*% ap- 
proaches zero. For large negative values of x it is 
assumed that F(x) approaches zero in the same way. 
Such initial distributions are said to be of the exponen- 
tial type. This type includes the exponential, normal, 
and log normal distributions. 
the probabilities ®,(x) and ;®,(x) converge with in- 


Under these conditions, 


creasing ” toward 
(x) = exp(—e~”), P(x) = 1 — exp(—e) (1) 
where 
y = a(x — u) (1’) 
is a reduced extreme value. The expressions a and u 
are scale and location parameters, respectively. The 
second derivatives of the probability functions vanish at 


x = u. Consequently the parameter wu is the most 


probable extreme value. The mean and _ standard 
deviation of the extremes are 
r= uF y/a, s=n7/(vV6a) (4 0.3¢¢22 (2) 


where the upper (lower) sign applies to the maxima 
(minima). Let .V be the number of observed extremes. 


Then the two parameters may be obtained from the 
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TABLE | 
Expected Means and Standard Deviations of Reduced Extremes 





N oN 

20 0 1.0628 
21 0 1.0696 
22 0.! 1.0754 
23 0.5283 1.0811 
24 0.5296 1.0864 
25 0.5309 1.0915 
26 0.5320 1.0961 
27 0.5332 1.1004 
28 0.5343 1.1047 
29 0.5353 1. 1086 
30 0.5362 1.1124 
31 0.5371 1.1159 
32 0.5380 1.1193 
33 0.5388 1.1226 
34 0.5396 1.12355 
35 0.5408 1.1285 
36 0.5410 1.1313 
37 0.5418 1.1339 
38 0.5424 1. 1363 
39 0.5430 1.1388 
40) 0.5436 1.14138 
{ 0.5442 1.1436 
42 0.5448 1.1458 
43 0.5453 1.1480 
44 0.5458 1.1499 
45 0.5463 1.1519 
16 0.5468 1.1538 
47 0.5473 1.1557 
48 0.5477 1.1574 
49 0.5481 1.1590 
50 0.5485 1. 1607 
51 0.5489 1.1623 
52 0.5498 1.1638 
53 0.5497 1.1653 
54 0.5501 1. 1667 


extreme observations alone without any knowledge of 
the initial distribution by the modified method of 
moments with the help of the equations 


(l/a) = (s/oy), u= & F W/a (33) 

The difference in sign arises because the asymptotic 
distribution of the reduced smallest value is the mirror 
image of the asymptotic distribution of the reduced 
largest value. The numerical values of ty and ox, 
called the expected extremes and standard deviations, 
are given in Table | as functions of .V. 

Since the reduction of Eq. (2) is linear, a probability 
paper for extreme values can be constructed and used 
in the manner described in previous publications.* ° 
First, the .V largest observations are ordered in increas- 
ing magnitude, and a rank m is ascribed to each ob- 
This procedure is used in order to conserve 
Grouping should 


servation. 
as much information as is available. 
be used only if the number of observations is very large. 
If the same observation happens several, say /, times, 
the geometric mean of the first rank m and the last rank 


m+p—1 
/ 
m= Vm(m + p — 1) (4) 


is attributed to the p observed identical values. This 
rule is not applied if the last observation occurs p times. 
In this case both the (V — p + 1)th and the Nth 
observation are employed to preserve the number of 


observations. The plotting positions are obtained by 





V Van oN 

5d 0.5504 1.1681 
56 0.5508 1. 1696 
57 0.5511 1.1708 
58 0.5515 L. E23 
59 0.5518 1.1734 
60 0.5521 1.1747 
62 0). 5527 1.1770 
64 0.5533 1.1793 
66 0.5538 1.1814 
68 0.5543 1.1834 
70 0.5548 1.1854 
72 0.5552 1.1873 
74 0.5557 1.1890 
76 0.5561 1.1906 
78 0.5565 1.1923 
80 0.5569 1.1938 
82 0.5572 1.1953 
84 0.5576 1. 1967 
86 0.5580 1.1980 
88 0.5583 1.1994 
90 0. 5586 1.2007 
92 0.5589 1.2020 
94 0.5592 1.2032 
96 0.5595 1.2044 
98 0 5598 1.2055 
100 0.5600 1.2065 
150 0.5646 1.2253 
200 0.5672 1.2360 
250 0.5688 1.2429 
300 0.5699 1.2479 
400 0.5714 1.2545 
500 0.5724 1.2588 
750 0.5738 1.2651 
1000 0.5745 1.2685 
x 0.5772 1.2826 


dividing the ranks (or their geometric means) by V + 1. 
If the theory holds, the observed values plotted on 
extreme probability paper should be scattered about 
the line of Eq. (1) which is conveniently written 


x=u+ty/a (9 


where the parameters “ and 1/a@ are obtained from 
Eq. (3). 

The calculation of the parameters may be avoided by 
the use of an alternative general procedure valid for 
large samples. Instead of the extremes x consider the 
standardized extremes: 


1, = (x = 2)/s (6) 


or 
t. = 0.77970y ¥ 0.45005 (7) 


For large samples all standardized extremes which 
follow the probability law of Eq. (1) should lie on a 
single line for ¢, as a function of y, with common mode 
0.45005, and the same slope = 0.77970. For small 
samples this technique may be improved by using Jy 
instead of y and 1/oy instead of +/6/r. Then Eq. (7) 
becomes 


-—/ 


i = (y a Vw) ‘ON (4) 
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and the corresponding line for ¢, as a function of y holds 
for all samples of size VV which follow the probability 
law of Eq. (1) whatever the extremes may be. 

The extremal paper shows the largest (or smallest 
values x on the ordinate and the probability (x) |or 
$(x) | and the reduced extreme value y on the abscissa. 
In addition an upper horizontal scale shows the so- 
called return periods 7(x). They are obtained in the 
following way: The probability of a value equal to or 
larger than xis 1 — ®(x). Its reciprocal is the expected 
number of observations necessary to obtain such a 
value once. If the observations (as is frequently the 
case With extremes) are made in regular intervals of 
time, 


Tix) = 1/]1 — ®(x) (S) 


is again a time and therefore called the return period. 
Since the return periods 7 are plotted on a parallel to 
the horizontal y axis, the extreme values x plotted 
vertically can be read off as functions of the return 
period. 

For large values of x, Eq. (5) becomes from Eqs. (8) 
and (1) 


x=utilgT/a (9) 
Consequently, 1/a@ becomes the rate of increase of the 
variate considered as a function of the logarithm of the 
return period. Once the line has been fitted, it can be 
extrapolated graphically, and the most probable ex- 
treme values corresponding to given return periods 7 
can be obtained. 

An important use of this theory in aeronautics is the 
extrapolation of the line for a certain number of flight 
Each record corre- 
There- 


hours instead of return periods. 
sponds to a certain number + of hours of flight. 
fore the return periods are multiplied by this number, 
starting with the return period of the mode u which is 
from Eqs. (8) and (1) 


T(u) 1) (10) 


w= 
= 1,582 


This leads to the consideration of the variates as func- 
tions of the hours of flight. If, instead, mileage is used 
as scale, the products 7-7’ must again be multiplied by 
an average speed. 

In many cases the scatter of the observations about 
the theoretical line Eq. (5) is so fine that the theory can 
be accepted on the basis of a visual inspection alone. 


TABLE 2 
Control Curves Near the Median 


V O(x)|1 — P(x) 

y ’'(x 
—0.5 1.2438 

0.0 1.3108 
+0.5 1.5057 

1.0 1.8126 

1.5 2.2408 

2.6 2.8129 
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In other cases, a doubt may exist whether a certain 
deviation is still compatible with the theory. To de 
cide the question of how far the observations can be 
the 


curves are constructed showing upper and lower limits 


scattered without invalidating theory, control 
within which the values can vary with a prescribed 
probability of, say, 0.67. This level is chosen because 
it corresponds to the one sigma limit for the normal 
For the the control 


curves,’ consider first the values x,, in the neighborhood 


distribution. construction of 
of the median, i.e., those values which are located in an 
interval, say, 0.15 << ®@ < 0.85. <A general theorem of 
statistics’ affirms that such values are normally dis- 
tributed about their mean, which is the value located 
on the straight line, and have standard errors 


I V P(x)[1 — P(x) ] 
aVN ?’(x) 


The expression in parentheses can be calculated as a 


function of the reduced variate y. This is done in 


Table 2. 
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The control curves are obtained by multiplication of 
the values given in Table 2 by the factor 1/(av.V) 
which is obtained from Eq. (3). The values s, are 
added to, or subtracted from, the theoretical values x 
located on the straight line at the corresponding values 
of y. 

Since oy (Table 
control curves can be intuitively interpreted as lines 
nearly parallel to the theoretical line of Eq. (5) at dis- 
tances which are of the order of the standard error of 
Obviously the spread of the control curves 


varies only from 1.0 to 1.3, the 


the mean. 
shrinks with increasing sample size NV, and increases 
with increasing standard deviations. 
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Eq. (11) cannot be applied to the extremes of the ex. 
tremes, since these mth values are distributed about 
their means not as normal values but as extremes. 
This leads to the control intervals 


1 = (0.754 a}, 
Axy—-2 = (0.589/a) (12 


Axy = (1.141/a), Axy 


valid for the largest and second and third extreme 
values. As previously, the A factors are added to, and 
subtracted from, the theoretical values on the straight 
line, and finally the two or three outstanding points are 
joined graphically to the previous control curves. If 
the number of extremes, .V, is small, it is sufficient to use 
Axy. If Nis so large that the largest value is far away 
from the control curves valid for the median, it is ad- 
visable to use Avy —,; and Axy—» also. Thus the control 
curves cover the whole observed range with the excep 
tion of the minima of the largest and the maxima of the 
smallest values which are not very interesting. 

To control the forecast, parallels to the theoretical 
line are drawn at a distance of 1.141/a. Since the re 
turn period scale as a function of y is approximately 
logarithmic, the control curves expand considerably for 
large values of x. This is not surprising; the longer the 
forecast wanted, the less reliable must it be. 

The same methods may also be used for smallest 
values if the variate is traced in decreasing order. This 
leads to an easy comparison of the increase of the 
largest and the decrease of the smallest value. 


(2) Extreme Winds 


Before applying these methods to aeronautical 
problems, it seems reasonable to ask whether such an 
unreliable feature as windspeed, and particularly its 
extreme, can be submitted to a statistical analysis. 
An affirmative answer has been given’ where the ex- 
treme windspeeds for the years 1912-1948, observed at 
Des Moines, Iowa, were studied by the theory of ex- 
treme values. 

A. Court! analyzed the strongest 5-minute windspeed 
(1912-1948) at 25 places in the United States and 
found good agreement with the theory of extreme 
values. The data of Table 3 taken from Court show 


the annual maxima of the greatest speed over a 5- 


TABLE 3 
Maximum Annual (5-Minute) Windspeed, Eureka ( 
1912-1948 


California 


M.P.H. Frequency Plotting Position 
30 po 0.046 
31 3 0.13 
32 1 0.18 
34 6 0.27 
35 10 0.47 
36 ] 0.63 
37 4 0.70 
38 4 0.802 
39 1 0.868 
40 l 0.895 
41 1 0.922 
46 2 0.947-0 .974 
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minute interval at a height of 10 feet. They are 
traced in Fig. 1 together with the theoretical line and 
the control curves. The estimates of the parameters 
and data for the construction of control curves are 
given in Tables 11 and 12. The fit is satisfactory. 
The single point outside the control curves may be due 
to the grouping of the observations. 

In the same way A.I. Johnson* studied maximum wind 
velocities and maximum mean wind velocities per hour 
at 25 Swedish stations and obtained a very good fit of 
the theory of extreme values. 

Next studied will be the influence of gusts on aircraft. 

An airplane flying through gusty air encounters large 
loads. In the design of airplanes large gust velocities 
are of particular interest since they are likely to cause 
structural damage. Consequently the statistical in- 
terpretation and forecast of gust velocities is an im- 
portant problem. 

Harry Press!! showed that the initial distribution 
from which effective gust velocities are derived is ex- 
ponential. Since the asymptotic distribution of the 
largest value is reached for this initial distribution even 
for small sample sizes, a good fit of the asymptotic 
theory to the observed extremes may be expected. 
Press applied this theory to observations derived from 
specific test transport 
operations, and proved that it leads to better results 
than the empirical approach, the fitting of a Pearson 
Since each maximum is taken 


conditions and commercial 


type III distribution. 
from a different traverse, the observations are inde- 
pendent, and the application of the theory of extreme 
values is proper. 

Press uses analytical methods requiring numerical 
calculations which can easily be avoided by the use of 
extreme probability The maxima of the 
effective gust velocities, in ft. per sec., for 485 traverses of 
thunderstorms, obtained from the 1946 Operation of the 
U. S. Thunderstorm Project, are given in Table 4. 
Columns | t are plotted in Fig. 2 on extreme 
probability paper. The calculation of the parameters 
from Eq. (3) is shown in Table 11. The factors ¥45 and 
o43, are obtained by linear interpolation from Table 1. 
The parameters u and 1/a calculated here differ slightly 


paper. 


and 


S 
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TABLE 4 
Effective Gust Velocities 
Variate, Frequency, Mean Rank, Plotting Position, 
x p ml m/486 
3 4 2 OOO 0.0041 
5 11 8.660 0.0178 
7 27 25.82 0.0531 
ie) 18 62.24 0.128 
11 62 117.6 0.242 
13 58 179.2 0.369 
15 55 236.5 0.487 
17 60 294 0 0 
19 61 354.7 0 
21 36 $04.1 0 
23 17 130.9 0 
25 18 148.4 0 
7 8 161.5 0 
29 rf 1690 ) 
31 6 475.5 0 
33 3 1800 0 
35 l 182.0 0.9918 
37 2 $83.5 0.9949 
39 l $85 0 0.9979 
TABLE 5 


Maximum Airspeeds 


Maximum 


Airspeed, Mean Plotting 
Vmax. Frequency, Rank, Rank, Position, 
m.p.h p m m m/118 

152 12 1-12 3.46 0.029 
157 18 13-30 19.75 0.167 
162 22 31-52 10.15 0.340 
167 23 53-75 63.05 0.534 
112 17 76-92 82 62 0.709 
117 16 93-108 100,22 0.849 
182 6 109-114 111.47 0.945 
192 ] 115 115 0.975 
197 l 116 116 0.983 
197 ] 117 117 0.9915 


from the values given by Press. The fit of the theo- 


retical line 


x = 12.807 + 4.920y 


to the observations is so good that only the ends of the 
control curves given by Table 12 are traced in Fig. 2. 

For maximum airspeeds, J’,,,,, two series, given by 
W. G. Walker and Roy Steiner,'” 
They are composed of 8 sets of observations comprising 


will be considered. 


$11 records and representing 65,210 flight record hours 


TABLE 6 
Maximum Acceleration Increments 


Plotting 


Variate, Rank, Position, 
Anmax m m/27 
0.50 l 0.037 
0.61 2 0.074 
0.65 3 0.111 
0.66 } 0.148 
0.70 5 0.185 
0.73 6 0.23992 
0.80 7 0.259 
0.80 8 0.296 
0.82 9g 0.333 
0.91 10 0.370 
0.91 11 0.407 
0.98 12 0.444 
1.00 13 0.481 


Plotting 


Variate, Rank, Position, 
An max m m 27 
1.00 14 0.518 
1.00 15 0.536 
1.05 16 0.593 
1.05 17 0.630 
1.10 18 0.667 
1.20 19 0.704 
1.25 20 0.741 
1.30 21 0.778 
1.30 22 0.815 
1.35 23 0.852 
1.40 24 0.889 
1.65 25 0.9259 
1.80 26 0.9630 
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in commercial transport aircraft Operations in various 
parts of the world during 1933-1941. The data in 
lable 5 are taken from column C-III of the period I 
operations. 

Estimates of the parameters and data for the control 
curves are shown in Tables 11 and 12. The observa- 
tions are plotted along with the theoretical line in Fig. 3. 
It is seen that the fit is a good one. 

Large samples of observations on IJ’,,,,, taken under 
essentially the same conditions seem to be very scarce. 
In order to get a large sample, the standardizing tech- 
7) has been applied to the individual sets 
The re- 


nique of Eq. ( 
of Vax data for the period I observations. 
sulting sets of standardized observations have been 
pooled providing a set of 411 observations. The esti- 
mates of the parameters and data for the control curves 
are shown in Tables 11 and 12. The observations and 
the theoretical line plotted in Fig. + show a good fit. 
Since the data are pooled the use of a common mean and 
standard deviation for the purpose of extrapolation is 
not realistic. Therefore no scale for V2, is intro- 
duced. The conditions under which pooled extremes 
follow the double exponential law of Eq. (1) are un- 
known. 

Figs. 3 and 
maximum air speeds follow the theory of extreme 


are favorable to the hypothesis that 


values. However, an analysis of S82 records from 
Dayton did not lead to a satisfactory fit. Therefore 
this application requires further study. 


(3) Extreme Gust Accelerations 


H. Press!! gives 26 records for maximum gust acceler- 
measured in g units, each 
250 hours of flight. The 


ation increments An,,,,, 
record representing roughly 
data are given in Table 6. The observations are 
plotted in Fig. 5 on extreme probability paper. The 
details of the calculation of the straight line and of the 
control curves are given in Tables 11 and 12. The 
parameters obtained differ slightly from those calcu- 


lated by Press. The theoretical line 
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TABLE 7 
Extreme Acceleration Increments 
Plotting 
Variate, Rank, Position, 
Anmax m m/l 
0.65 l 0.032 
0.75 2 0.064 
0.85 3-4 0.112 
0.95 5-6 0.177 
1.05 7-10 0.270 
1.15 11-14 0.400 
1.25 15-20 0.559 
1.35 21-24 0.72 
1.45 25-26 0.823 
1.75 27-29 0.9038 
1.95 >) 0 O677 
x = 0.865 + 0.292y 


leads to an excellent fit. 


Press compared the total miles of flight with the 


estimated number of miles required to exceed the largest 
value of acceleration increment actually measured at 
each test. The agreement between the observed and 
theoretical values is satisfactory for the extreme value 
theory, and unsatisfactory for the Pearson type III 
distribution previously used to analyze these data. 
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TABLE 8 
Extreme Accelerations Am max 


Acceleration Airline A-—Prewar Airline A—-War 





Increments Plotting Plotting 
in g, Rank, Position, Rank, Position, 
Anmax m m/ol m m/41 
0.35 
0.45 
0.55 
0.65 l 0.032 
0.75 2 0.065 1-2 0.034 
0.85 3-4 0.112 eben 
0.95 5-6 0.177 3-4 0.084 

1.05 7-10 0.270 5-9 0.164 
3 11-14 0.400 10-11 0.256 
1.2 15-20 0.559 12-16 0.338 
1.35 21-24 0.724 17-20 0.450 
1.45 25-26 0.823 21-25 0.559 
1.55 26 0.634 
1.65 27-30 0.694 
1.75 27-29 0.903 31 0.756 
1.85 32-34 0.802 
1.95 30 0.968 35-37 0.877 
: 38 0.927 
2 35 39 0.951 
9 45 40) 0.9756 


Peiser and Wilkerson" analyze 30 observations of 
extreme acceleration increments by the Pearson dis- 
The data 


tribution. Their data are given in Table 7. 


are traced on extreme probability paper in Fig. 6. The 


calculations for the parameters and for the control 


curves are shown in Tables 11 and 12. The line 


x = 1.083 + 0.267y 


is traced in Fig. 6 together with the forecast obtained by 
Peiser and Wilkerson with the help of Pearson's type 
III distribution. This forecast falls short of the large 
observations and of the forecast obtained by the theory 
of extreme values. 

Table S shows four similar series given by Peiser.' 
They are traced in Fig. 7 together with the theoretical 


4 


lines obtained from the parameters given in Table 11. 
The fit is so good that there is no need to calculate the 
control curves. Of course, the best agreement between 
theory and observations is reached for Airline B since 
N = 74 is quite a large sample. 

In these observations the largest accelerations sur- 
passing | g and the smallest accelerations falling short 
of 1 g are lumped together. The accelerations | — A 
below unity are counted as if they were above unity 
by the same amount A. For commercial transport 
operation the accelerations are predominantly due to 
gust, and the available data indicate that the distribu- 
tion of positive and negative (up and down) gust 
velocities and accelerations are symmetrical. Under 
these conditions the artificial lengthening of the record 
is permissible. 

However, in general gust and maneuver types of 
acceleration must be separated. The symmetrical 
character of accelerations does not apply to most mili- 
tary type airplanes, since the accelerations involved 
here are largely pilot-induced in rapid maneuvers and 


Airline B Airline C 
Plotting Plotting 
Rank, Position, Rank, Position, 
m m/75 m m (23 
l 0.048 
2-3 0.106 
l 0.013 1-5 0.194 
2-4 0.038 6 0 261 
5 0.067 7-9 0.342 
6-14 0.122 10-12 0.475 
15-20 0.231 13-14 0 585 
21-28 0.322 15-17 0. 695 
29-42 0.465 18-19 0 S804 
43-51 0.624 
52-57 0.721 
58-62 0.795 20-22 0.870, 0.9565 
63-65 0.851 
66-67 0.886 
68-69 0.915 
70-71 0.940 
72 0.9600 
73-74 0.9733, 0.9866 


are characterized by large and more frequent positive 
than negative accelerations. The maxima and minima 
of accelerations observed in Dayton and given in 
Tables 9 and 10 are suggestive of maneuver rather than 
of gust accelerations. Since the distributions of the 
accelerations below and above unity differ and since the 
summation of two different distributions of largest 
values does not lead to another distribution of largest 


TABLE 9 
Maxima of Acceleration (>1.2 g, Dayton 


Plotting 


Variate, Frequency, Rank Mean Rank, Position, 


x p m m m/78 
1.20 4 1-4 2 0.026 
1.23 6 5-10 1.7 0.091 
1.30 10 11-20 14.85 0.191 
1.35 7 21-27 23.81 0.305 
1.40 10 28-37 31.19 0.400 
1.45 11 38-48 B.3 0.548 
1.50 10 49-58 52.32 0.671 
1.55 | 59 59 0.756 
1.60 ~ 60-67 63.4 0.812 
1.65 4 68-71 69.5 0.891 
1.70 2 72-73 72.5 0.9295 
1.75 2 74-75 74.5 0.9551 
1.400 l 76 76 0.9744 
20) l 77 ej 0). 9872 

TABLE 10 


Minima of Acceleration (<0.8 g, Dayton 


Plotting 


Variate, Frequency, Rank, Position, 
v p m m /66 
0.80 7 1-7 0.040 
0.75 6 8-13 0.155 
0.70 16 14-29 0.305 
0.65 11 30-40 0.525 
0.60 8 41-48 0.672 
0.55 8 49-56 0.794 
0.50 6 57-62 0.901 
0.45 2 63-64 0.9621 
0.35 l 65 0.9848 
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values, the data cannot be combined, and the two ex- 
treme accelerations have to be studied separately. 
The calculations for the parameters and for the control 
curves for both distributions are shown in Tables 11 and 
12. The observed frequencies are traced in Fig. § 
together with the theoretical straight lines 

x = 1.370 + 0.139y 
for the maxima, and 
0.692 — 0.085y 


x= 


for the minima of the accelerations. 
To make the two series comparable, the minima of 
the accelerations are traced in descending order such 
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the accelerations from 


larger than the mode of the absolute smallest deviation 


largest deviation of unity jis 
from unit’, and the spread of the first distribution js 
larger than the spread of the second. Fig. 8 shows that 
the largest accelerations increase more quickly with the 
return period, and consequently with the hours of flight, 
than the smallest accelerations decrease. 
CONCLUSIONS 

A study of extreme values in aeronautics is greatly 
handicapped by lack of data. A large sample of ex- 
tremes taken under essentially the same conditions is 
almost unheard of. The use of the asymptotic theory 
requires the existence of a fairly large sample of V 





that 0.8 corresponds to 1.2, etc. The mode of the extremes taken from subsamples of large size n, thus 
TABLE 11 
Estimation of the Parameters 
Table 3 i 5 6 7 
Effective Maximum Extreme 
Extreme Gust Acceleration Acceleration 
Variate Windspeed Velocity Maximum Airspeeds Increment Increments 
Dimension m.p.h f.p.s m.p.h m.p.h g @ 
Source A. Court NACA 1926 NACA 2625 NACA 2625 NACA 1926 NACA 807 
Number of observations N 38 185 117 111 of 20) 
Mean t 35.595 15.623 166.5 0 1.9200 1 2267 
Standard deviation s 3.7513 6.1899 9.66 1 0.3202 () 2074 
Factor from Table 1 oN 1. 1363 1.2582 2.13 1.2549 1.0961 1.1124 
Parameter l/a 3.277 4.9196 7.966 0.7969 0.292] (). 2674 
Factor from Table 1 VN 0.5424 0.5723 0.5616 0.5716 0.53820 0.5362 
Mode u 33.819 12.807 162.1 —(. 4554 0.86559 1.0833 
Factor for control curves 1/(aV NV) 0.539 0.223 0.736 0.0392 0.0573 0.0488 
Figure 1 2 3 } 5 6 
TABLE 11 (Continued) 
Fstimation of the Parameters 
Table 8 8 8 8 9 10) 
Extreme Extreme Extreme Extreme 
Acceleration Acceleration Acceleration Acceleration Largest Smallest 
Variate Increment Increment Increment Increment Accelerations Accelerations 
Dimension g g g g g g 
Source NACA 1142 NACA 1142 NACA 1142 NACA 1142 Dayton Dayton 
Number of observations A 30 1) 74 22 77 65 
Mean t 1.2267 1.4500 1.2054 0. S864 1.4474 0.6446 
Standard deviation s 0.2974 0.3988 0.3378 0.3156 0.1627 0.1008 
Factor from Table 1 oN 1.1124 1.1418 1.1890 1.0654 1.1692 1.1824 
Parameter l/a 0). 2674 0.3494 (). 2841 0.2962 0.1391 0) O854 
Factor from Table | Vy 0.5362 0.5486 0. 5557 0.5268 0.5563 0.5538 
Mode u 1.0833 1.2601 1.0475 0.7303 1.3700 0.6918 
Factor for control curves 1/(aV NV) 0.0467 0.0551 0.9321 0.0632 0.0159 0.0106 
Figure 4 7 7 7 8 8 
TABLE 12 
Control Curve Calculations 
Reduced Gust Accelerations Accelerations 
Variate Windspeed Velocity Maximum Airspeed Maximum Extreme Maximum Minimum 
NACA NACA NACA NACA NACA 
y Court 1926 2625 2625 1926 807 Dayton Dayton 
—0.5 0.67 0.28 0.92 0.049 0.071 0.061 0.020 0.0135 
0 0.7] 0.29 0.96 0.051 0.075 0.064 0.021 0.014 
+0.5 0.81 0.3 L.# 0.059 0.086 0.073 0.024 0.016 
1.0 0.98 0.40 1.33 0.071 0.104 0.088 0.029 0.019 
1.5 1.2 0.50 1.65 0.088 0.128 0.109 0.036 0.024 
2.0 1.52 0.63 2.07 0.111 0.161 0.137 0.045 0.030 
AXn ; 2.90 +. 70 0.47 0.082 0.050 
AXn-1 2.47 3.71 6.00 0.600 0.202 0.105 0.064 
AXn 3.74 5.61 9.1 0.91 0.333 0.305 0.159 0.097 
Table 3 4 5 6 7 9 10 
Figure 1 r 3 4 5 6 8 8 
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EXTREME VALUES 


requiring very large numbers of observations. Some 


for the most part, 


researchers in aeronautics have, 
taken data with the idea of performing a classical 
statistical analysis, and hence have had no particular 
interest in recording extremes. Others have tried to 
analyze maximum gust velocities, for example, and 
have been forced to use pooled observations on several 
types of aircraft under as many different operating 
Still others have collected extreme value 


without 


conditions. 
data under various conditions noting the 
means and standard deviations of the observed sam- 
ples. This shortcoming has not always been avoidable, 
as some of the recording instruments used record, by 
design, only an envelope of extremes. 

It is hoped that future plans for collecting extreme 
values data will be designed to obtain observations 
under the same operating conditions. Where it is not 
possible to obtain such data, some effort should be 
made to estimate the mean and standard deviation of 
the underlying population from which each sample is 
drawn. This will permit intelligent adjustment of the 
data if pooling is necessary. 

The heuristic method of analyzing maximum air 
speeds and largest accelerations by fitting a Pearson 
type III distribution does not take into account the 
most important intrinsic property that the observed 
are extremes. The method used requires the 
third moment which is strongly 


values 
calculation of the 
aifected by errors of random sampling, and which be- 
comes reliable only for sample sizes which are not 
available. 

The method outlined here is logically adequate for 
values, avoids the calculation of the third 
is easier to handle by the routine procedure 
shown and better fit to the 
Graphical techniques lead to a criterion for the good- 


extreme 
moment, 
gives a observations. 
ness of fit and to a forecast. 
For the military type of 
lengthening of records by combining extreme negative 


airplanes the artificial 


and extreme positive accelerations leads to a serious 
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underestimation of the dangers involved in the occur- 
rence of extreme values. In general the new theory 
leads to a larger forecast of extremes. 

In the cases studied here, the fit of the theory of ex 
treme values is very good. Therefore it is advocated 
that the asymptotic theory of extreme values be used 
for the analysis of gust velocities, maximum airspeeds 
and largest and smallest accelerations. In particular 
this method should be used to replace the empirical 
procedures used up to now in the construction of the 


V-G diagram. 
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The Effect of Leading-Edge Bluntness on a Laminar 
Supersonic Boundary Laver 
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cone show little 
effect of tip blunting on the laminar boundary layer. 


Comparable experiments on a 15° 
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Charts for the Calculation of the Critical 
Compressive Stress for Local Instability 
of Columns with Hat Sections 


CHRISTIAN J. VAN DER Maast 


Fairchild Engine and Airplane Corporation 


SUMMARY 


Charts are presented for the calculation of the critical com 
pressive stress——the stress at which local instability occurs—for 
olumns with hat section, and an outline for the use of the charts 
sgiven. This method can be used to calculate the critical com 
pressive stress above, as well as within the elastic range 


INTRODUCTION 


ee MAXIMUM aerodynamic efficiency of a wing, the 
compression panels in the wing should not buckle, 


even when the applied load reaches high values. This 
requirement makes it necessary for aircraft designers to 


determine the critical compressive stress at which local 
instability occurs. Local instability of a column is de- 
fined as any type of instability in which the cross sec- 
tions are distorted in their own planes but are not 
translated or rotated. 

As an approach to solving these problems for columns 
with hat section, a method of calculating the elastic 
theoretical buckling stress for local instability of hat sec- 
tions is presented. This method is extended to the cal- 
culation of the critical compressive stress above the 
elastic range through the introduction of the plasticity 
correction factor 7. 

This report includes charts of the buckling co- 
efficient, along with tables of the values used in prepar- 
ing the charts. An example of the use of the charts is 
also presented. 


SYMBOLS 


= width of plate element of section 


= thickness of plate element of section 


E = Young’s modulus of elasticity 
Etan = tangent modulus of elasticity 
Esee == secand modulus of elasticity 

m = Poisson's ratio (taken as 0.3) 
er = critical compressive stress 

a = compressive yield stress 


= plasticity correction factor 
k = nondimensional buckling stress coefficient dependent 


upon relative dimensions of cross section 


Received November 17, 1953 
*The author wishes to thank L. T. Waterman for his help 


ful cooperation in setting up this problem on the I.B.M. Card 


Programmed Calculator 
+ Stress Technologist, Fairchild Aircraft Division. 


399 


K = () 903k 

o = stiffness factor 

& carry-over factor 

S, C = far end fixed 

sl Cc! far end elastically restrained 

sul far end free 

SY far end supported and subjected to sinusoidally dis 
tributed moment equal and opposite to moment 
applied at near end 

S = S/(D/b 


Suser pts 


F = flange 
W = web 
t = top 
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_W 1.0 - COMPRESS ION BUCKLING CONSTANT 
‘ = VK; FOR HAT SECTIONS 
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Compression buckling constant +~/ A; for hat sections, 
continued. 
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Compression buckling constant ~ A; for hat sections, Fic. 5. 
eontinued. 


Fic. 4. 
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1.6 T T T 


tE 10, tWe.st 


is shown in Fig. 1. The critical compressive stress may 1.2 
be computed from the formula } si 
~ i 
xs rE i \* 
= - ( : ) = KOE ) (1 es a 
n "1201 + bu”) \b, b, Be es Sag ee |bp/bw 


where Sg ie, i 


/ Dy n = 0.5 + 0.5V 0.25 + 0.75(Eyan/Esec) (2) 4 Te 
‘ ie, 


K, = 0.9032, 


CALCULATION OF CRITICAL COMPRESSIVE STRESS 


The idealized cross section considered in this report 


1 








For a given section all the quantities on the right-hand 





ie : 1 4 1 

side of Eq. (1) are known, except the value of the co- 0 10 12 1.4 1.6 
efficient k,, This unknown value may be determined by/by 

from the appropriate chart (Figs. 3-6) after the neces- 


o 


. . : : Fic. 6. Compression buckling const: A, for hat s s 
sary dimension ratios have been determined. a Se ree ilies: 


DISCUSSION OF CHARTS The series criterion for stability as presented in ref 


; ‘ ‘ : erence 2, was selected 
Values of &, used in the preparation of the charts were 


























computed by an application of the principles of moment r = (Sw*Cw*)/[(Sw + Sp'")(Sw + S'%)] = 1 (3) 
distribution to the stability of thin plates. 
Because of the symmetry of the section about the ” 
vertical axis, only two joints may be considered for r = (Sy?Cyw?)/} [Sw + TA x 
stability consideration. [Sw t+ (5: cal} = 1 
Table 1 
Sample Calculation of kt 
kt At/ot |ow/de | bp/ow | tw/te |tr/ty [5,27 | rw/ow | ky | 5S [S%c? |apfoe | kp |SpUT] d, Dp r ke | (ky) r= 
[cove # 105 104 | 100 | 102 | 1@ | 103 | 213 | 106 | lo7 | 21% |215 | 108 | 109] 220 | 120 | a | 222 | 105 
1.020 | .900 |.800 |.600 | .510 |1.000 |1.225 | 1.125 ]2.510] 1.193 | .355 |1.875 |.904|-.631 | 8.582 | .1412] 293 
1.040 1.228 2.559 | 1.184 | .363 921 | -.717 | 8.589 |-.0125 | -3.662 | 1.0332 
1.0331 1.227 2.542 | 1.187 | .360 .915 |-.687 |8.587 | .0427| .983 | 1.0332 | 1.0332 
1.020 | .950 1.180 | 1.1875]2.510| 1.162 | .356 |1.979 |.904 |-.666 |8.278 | .0506] .350 
12040 1.183 2.559 | 1.153 | .364 921 | -.748 |8.286 |-.0937 |- .468 | 1.0221 
1.0211 1.180 2.522 | 1.161 | .357 -904 | -.671 | 8.279 | .0432 | .9967 | 1.0212 | 1.0212 
1.020 | 1.000 1.132 | 1.250 |2.510] 1.136 |.355 | 2.083.904 | -.676 |7.966 | .0101 | 4.409 
1.040 1.128 2.559 | 1.120 | .362 .921 | -.751 |7.930 |-.1232 |- .370| 1.0147 
1.0147 1.134 2.497 | 1.138 | 353 899 | -.657 17.975 | .0437 | 1.0135 | 1.0146 
1.0146 1.134 2.497 | 1.138 | .353 .899 |-.656 |7.976 | .0443 | 1.0002 | 1.0146 | 1.0146 
1.020 | 1.100 1.004 | 1.375 |2.510] 1.100 |.349 | 2.292 |.904 |-.647 |7.158 | .a224 | 2.173 
N 1.040 1.000 2.559 | 1.094 | .354 -921 | -.709 |7.125 |-.0870 |- .572 | 1.0151 
1.0151 1.006 2.498 | 1.102 | .348 -899 | -.632 17.166 | .0481 | 1.007 | 1.0150] 1.0150 
1.020 | 1.200 -909 | 1.500 2.510] 1.076 | .341 | 2.500].904 |-.589 | 6.560 | .0951 | .546 
1.04¢ -905 2.559 | 1.071 | .345 .921 | -.638 |6.531 | .0071 |7.497 | 1.c@298 
1.0298 -907 2.534 | 1.073 | .343 .912 |-.613 16.546 | .0526 | 0966 | 1.a298 | 1.02908 

































































Plot (ky) r = 1 against \¢/b_. ‘This gives: ky = 1.013 for At/>y = 1.050 


kf is obtained by interpolation: r(kye-b) =a 
ry kt) - rq kt 
Where: be. De. eal Bon 

r2-Tl 

So Thats a= rjkt) - rjb 


ons, 
: This gives (rel): Ky = a+b 
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Table 2 


Calculated Minimum Vaiues of Vky For Hat Sections 
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The step by step solution is as follows: 


(1) Determine the dimension ratios by },, by dy, 


ly f, and te ly. 
























































tr/ty = 1.0 
we yar (2) Assume A, 0,. 
Vie (3) Calculate Ay Oy (A, 61) (b,/ by and ) 
2 <0 by = (% bint bw) (by, bp). 
(4) Assume a series of values of &,. 
2 | a of | 5 | -6 | | 1.0 : ; oe 
(5) Determine kw = k,(by/b,)?(t,/tw)?; and ky = 
= 1.0 
tw/te ( at (by by (ly ty a 
3 T 2.276 53 2-275 3 7am | 5am (6) Read S,'¥ as a function of \, 6, and ,. 
“4 i . = 2.26 ‘ 2.238 | 3.008 (7) Read Sy and Sy?Cy* as a function of Ay by 
045 ad - - - - 2.046 and ky. 
a - - - 2.216 | 2.213 | 2.175 | 1.855 (S) Read S,'"! as a function of \p/by and kp. 
55 m Z s i, 2 2.050 . (9) Determine D,, D2, and (D;)(D2) [Eqs. (4) and 
“6 - 2.191 | 2.190 |] 2.187 | 2.1275 | 1.895 | 1.552 \) J. ; 
> ),)(Ds) (Item 9) ; ’Cw? (Item 7 
. ee eee ee ee ee ee (10) Plot (D,)(D.) (Item 9) and Su Cu Item i 
against k, (Item 4+) and note the value of k, for which the 
8 2.116 | 2.115 | 2.104 | 2.045 | 1.834 | 1.430 | 1.165 : 
two curves intersect (7 = 1). 
a ae ae | eae ae | ee et (11) Repeat steps 2 through 10, assuming different 
1.0 1.944 | 1.933 1.876 | 1.699 1.478 | 1.146 933 values of rN b,. 
£4 1.695 | 1.678 | 1.602 | 1.425 | 1.234] .956 -778 (12) Plot values of k, for r = 1 (Item 10) against 
1.5 1.384 | 1.367 | 1.297 | 1.144 989 | .765 -623 \, &, to determine the minimum value of k,. This 
ty/ te = .79 
8 1.931 - 1.869 - 1.467 | 1.135 924 Table 2 (Cont'd) 
1.0 1.656 - 1.549 - 1.176 | .909 -740 
1.2 1.402 “ 1.301 “ -981 | .758 -617 A 
a bw 
164 ee - 1.119 - -841 | .650 +529 — 
be br/bw 
1.6 1.062 - 082 - -737 | .569 463 
2 | # | 6 | 8 | 1.0 
Where ty/te = +63 
S = S/(D/b) a = 2.075 = - = 
D = Et®/(12(1 — p?)] 24 - 2.041 - - - 
Z _ be t, \3 5 - 2.006 - ~ - 
(Seca - (.S;) ( )( - 
b; li .6 ” 1.945 = ~ = 
Write .8 | 1.712 | 1.575 | 1.277 | .ow | .736 
r= (Stl) [(D,) (D2) | = ] (4) a Pa 1.404 “- - - 
where 1.0 | 1.380 | 1.266 | .943 | .727 | .592 
o SIV . 
Dy = Sw + (St) rea (5a) 1.2 1.154 | 1.057 | .786 | .606 | .293 
: 6 II - 
D, = Sw + (Sr rea (5b) 1.4 .991 007 | .674 | .520 | .223 
Except for this criterion, the step by step procedure is 1.5 2 .8L7 = ‘i = 
very much as developed in reference 1: The criterion 
(3) must be found by trial and error for each db, while 1-6 Oe “795 he °455 _— 
the lowest load that satisfies Eq. (3) is the critical load tyy/t, = -51 
for the section. = 
rr . ° ’ i . . . 9 of o 
rhe appropriate values of S and C for any \/6 and k . sata ae 956 6d 
may be found in the tables of reference 4. Because of 1.0 1.140 1.037 -765 -590 - 480 
awkward dimension ratios, it is often necessary to inter- ° . os 638 . 400 
polate between tables of reference 4. Therefore, it was ” 7 a - ‘ 
decided to go back to the original formulas for S and C es 815 e741 547 422 0343 
as presented in reference 1 and to perform the calcula- 
1.6 713 | .6 | .479 | .369 | .300 


tions on the I.B.M. Card Programmed Calculator. 
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CHARTS FOR CALCULATION OF 
minimum value of &, is the point on the curve for the 
appropriate dimension ratios. 

|3) Repeat steps | through 12, assuming different 
dimension ratios. 

A sample calculation of k, is given in Table 1. Table 
2 gives the values of k, as calculated with the method 
mentioned above. They are plotted in Figs. 3-6. 

The values of k, were calculated on the assumption 
that there is no lateral movement of the flange (see Fig. 
2a). This assumption is not justified, however, when 
the flange is of such proportions that it has small flexural 
stiffness in the lateral direction. The general case of 
cross-sectional distortion is of the type shown in Fig. 
2b, and the value of k, for this case is less than the value 
calculated on the assumption of no lateral movement. 

When a hat section buckles, one of the elements of the 
cross section is primarily responsible for the instability, 
inasmuch as this element is no longer stable in itself and 
requires restraint from the other elements until the 
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section as a whole becomes unstable. The charts show 
which element of the cross section is being restrained 
against buckling by the other elements. The dashed 
lines on the charts show the points for which an element 
is starting to be destabilizing upon the section. 

It seems likely that the maximum efficiency for local 
instability is in the range of dimensions near the point at 
which buckling is equally due to two elements. 
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Application of Volterra Linear Integral Equations to the 
Numerical Solution of Vibration Problems—II 


(Continued from page 388) 


" 


P, = same as P,; except ];,” substituted for },’ 


O, = (h/2) (Yo’ + 2Y,’ +... + 2Y,-1’ + Y,’) 


Q. = (h/2) (Yo" + 2Y," +... + 2¥n” + Y,”) 
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Skin Friction and Heat Transfer 
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(1) SUMMARY 


This paper presents the derivation of compressible turbulent 
boundary layer formulas for skin friction and heat transfer co 
efficients for the case when mass is injected into or withdrawn 
from the boundary layer normal to the solid surface. Curves 
are given illustrating the effects of mass transfer on local and 
average skin friction coefficients, mean velocity distribution, and 
mean temperature distribution for various Mach Numbers, Reyn- 
olds Numbers, and wall temperature to free stream temperature 
ratios. A comparison is made with some recent incompressible 
turbulent boundary layer experimental results with mass trans- 


fer present, and good agreement is shown. 


(2) SYMBOLS AND NOMENCLATURE 


A? = (Te/Tw)|(y — 1)/2] Mo, see Eq. (19) 
B = (Tao/Tw) — 1+ A?, see Eq. (19) 
c = see Eq. (21) 
Cy = 2ry/pxol'*, local skin friction coefficient 
° 
. = x! f, c;dx, mean skin friction coefficient 
CH = local heat transfer coefficient, see Eq. (44) 
Cu = mean heat transfer coefficient, see Eq. (46) 
c = specific heat at constant pressure 
" = 2w/pexl’«W, local mass transfer coefficient 
C = x ‘- c,dx, mean mass transfer coefficient 
D = empirically determined constant, see Eq. (23) 
E = see Eq. (26) 
h = heat flow 
I = integral number n 
= coefficient of heat conductivity 
7 = laminar eddy coefficient of heat conductivity, see 
Eq. (4) 
K = mixture length proportionality constant 
l = mixture length,/ = Ay 
M = Mach Number 
p = static pressure 
q = k(d7T/dy), heat transfer 


R = perfect gas proportionality constant 
R, = pol woX/po, Reynolds Number 
absolute temperature 


u, 1 = x,y components of velocity, see Fig. | 

w = mass flow per unit area injected into the boundary 
laver at the solid surface, see Fig. 1 

x,y = rectangular coordinates, see Fig. 1 


g = u/l 
Greek Symbols 


ratio of specific heats 


7 = 
6 = boundary layer thickness 
6* = boundary layer displacement thickness, see Eq. (49) 
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q = turbulent (eddy) viscosity 

c = variable of integration 

6 = boundary layer momentum thickness, see Eq. (29 
Me = viscosity coefficient 

p = density 

T = shear stress 

w = exponent of viscosity-temperature relationship, se« 


Eq. (33) 


total or inherent 
w = wall value or value at y = 0 
value at edge of boundary layer at y = 6 


= free stream value 


(3) INTRODUCTION 


— AERODYNAMICS includes many problems 
dealing with the flow of a high velocity, high 
stagnation enthalpy gas over component surfaces under 
turbulent boundary layer conditions. The flow of hot 
exhaust gases through rocket nozzles, of heated air 
through supersonic-hypersonic wind tunnels, and of 
atmospheric air along the surface of high velocity air- 
craft and missiles include but a few examples of flows of 
this type. Adiabatic wall temperatures under these 
conditions can exceed temperature limitations of struc- 
tural materials commonly used in the manufacture of 
such devices. A possible method of reducing the con- 
vective heat transfer from the hot gases to the cool 
surface is that of forcing mass through the deliberately 
porous wall to diminish the convective heat transfer to 
the surface. This paper deals with the effect of a small 
amount of mass transfer on the heat transfer from a 
high-speed, turbulent, compressible boundary layer to 
a smooth, porous surface. 

Yuan,! Morduchow,? and others* have considered the 
corresponding problem for laminar boundary layers. 
This paper is intended to provide relations for use when 
turbulent boundary layer conditions might be expected 
to prevail. Since the mechanics of the compressible 
turbulent boundary layer with no mass transfer is at 
present inadequately described, the present analysis 
wants experimental evidence to establish its limitations. 
skin friction co- 
upon Reynolds 
to free 


Formulas for local and average 


dependent 
temperature 


efficients are presented 

Number, Mach Number, 
stream static temperature ratio, and the ratio of rate of 
rate. 


wall 


mass transfer to local free stream mass flow 
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FIG Boundary-layer flow along a plane surface with mass 
injection normal to the wall 


Figures illustrating the effect of mass transfer on these 
coefficients and the turbulent boundary layer velocity 
and temperature distributions are presented along with 
comparisons with recent experimental results. 


(4) THEORETICAL DEVELOPMEN1 


(a) Equations of Motion 


The equations of motion for the thin compressible 
turbulent boundary layer over a flat plate have been 
These equations apply without 


given by Van Driest.* 
However, the bound- 


modification in the present case. 
ary conditions for the case of mass transfer normal to 
the surface of the plate differ from those treated pre- 
viously by Van Driest and therein lies the difference be- 
tween his development and the present case. The basic 
equations of motion will be repeated here along with 
their boundary conditions in order to illustrate the ap- 
proximate solution developed in this paper. 

The x momentum equation is, including the viscous 
shear stress significant in the sublayer |Van Driest’s 
Eq. (16) including the term [u(07, Oy) J. 

[pu(On/Ov)|] + [pv (0% Oy)] = O7/Oy (1) 
where 
r = [u(Ou/Ov)| + [e(Ou/Ov)] = 
[u(Om, Ov)] — (pv)’u’ (2) 
The y momentum equation is 
Op oy = 0 (3) 
The energy equation is, including the viscous dissipa- 


tion term and the laminar heat conduction term 


}Van Driest’s Eq. (26) including terms u(0#/Oy)* and 
(0/ Oy) [k,(OT/ Oy) ]}. 


0(C,T) 4 0(C,T) 
u v = 
‘ Ox . oy 
° ( or, , OF ed on \* (4) 
) d 
dy \" oy dy wide 
where 
oT 
k = —C, (pv)’T” (5) 
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The continuity equation is 
O( pit ) O( pi 
p 4 f 
Ox Ov 


or, for the case of mass transfer at the wall 


O( pu . a 
pi Ww “ ay (4 
S00 On 


where w = mass flux per unit area at the surface in y 
direction. 
The boundary conditions for the present problem are 


(I) at the surface, Fig. 1 


y 0; T 1 (Sa 
ul 0 Sb 
p p Sc 
pi Z Sd 


(II) at the edge of the boundary laver 


\ 0 T 1 Ya 
ul l (9b 
p p (Ye 
pv p (9d 


Because of Eq. (5 
p = p = p (10 


It is immediately apparent that for the case of con 
stant specific heat, C,, and Prandtl Numbers (C,u)/k 
and (C,¢)/k equal to | that the particular integral to 
Eq. (4) relating T to #@ given by Crocco® and Van 
Driest‘ applies in the present case since the introduction 
of boundary condition (Sd) does apply in finding the 


integral to Eq. (4). This integral is [Van Driest’s Eq. 


(34) ] 
T - 


Fa Pealli- PB) 4 (jue) 
(#) an 


However, it will be shown in the following discussion 
that the solution to Eq. (1) will differ from that pro- 
posed by Van Driest due to the additional boundary 
condition (Sd) contributed by the mass transfer at the 


surface. 


(b) The Velocity Distribution 
An approximate velocity distribution (x,y) will be 
derived in this section in a manner analogous to that 
used previously by Prandtl,’ and von Karman’ for in- 
compressible flow, and Van Driest* for compressible 
flow—both without mass transfer. 
Consider Eq. (1) evaluated at the plate surface, y = 


0, and with pv, at y = 0 given by Eq. (Sd) 
w(O0u/Oy)y-0 = (Or/Oy), 


or 
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[(0/Oy)(7 — wit) |,-0 = O (12) 


Now, for small values of y at station x = £, expand 


about (£,0) in a Taylor's series, using Eq. (Sd) to obtain 


r(£,0) + 
[(0/Oyv)(r — wi) |y=ov + 


r(é&,v) — wuléy) = 


or, because of Eq. (12), near y = 0 


(13) 


T— WU = Ty 


Now Prandtl,’ voi Karman,® and Van Driest! were 
dealing with the case of no mass transfer when w = 0. 


Thus, they equated the shear stress 7 (x,y) to the wall 
stress, sought a relation between the derivatives of 7 
and 7 (x,y) and arrived at a velocity distribution law 
relating %, p, € ¥, Tw, Px and empirically determined 
An analogous procedure will 
[While it is 
Length 


constants of integration. 
be followed in the present discussion. 
generally agreed that the Prandtl Mixture 
Theory’ and the assumption implied in using Eq. (13) 
do not describe accurately the kinematics of the turbu- 
lent boundary layer, nevertheless collectively they are 
quite successful in yielding useful skin friction laws. 
See, for example, reference 10. |] 
From Eggs. (2) and (13) 


Tr = — (pv)'u’ + p(Ou/Oy) — wi 


Or, expanding the correlation 


Try = — pur’ —G p'u’ — p'u'v’ + 


u(Ou/Oy) — wih (14) 


Assume, as is consistent with incompressible turbulent 
boundary layer measurements, that within the thin 
turbulent boundary layer the terms 


—vpu'; — p'u'v': u(On Oy) 


are negligible compared to the remaining terms in Eq. 


(14). Then make use of the Prandtl mixing length, /, 
where 
= , aul” 
uy = —f 
dy 


With / = Ay Eq. (14) becomes 
Ty = pK*y?(du/dy)* — wit (15) 


Eq. (15) can be rearranged to obtain the differential 
equation for the mean velocity component “(y) 


dz r 
tn = dy/Ky (16) 
px 
| — = a - “| 
B 2 
where 
¢ = a/U.: ¢, = 27,/(p.U.*): ——— 


2w/(p.U.) 


(17) 
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and, from Eqs. (17) and the Equations of State, ) = 
pRT 
p./p = T/T. 
Calling 
A? = (T./T.)[(y — 1)/2)]M.?; B = (T../T,,) - 
1+ (T./Tw)[(y — 1)/2)M.? (19 


Eqs. (11), (18), and (19) can be used with Eq. (16) to 


obtain 
aero 
dz AS =) (1 + Cs)\(1 + Bs — A2z?)]'* = 
dy/Ky (20 
where 


Integrate Eq. (20) to obtain 


1/K log, y = 


» / T. Ys ; 
S dz} T | (1 + Cs)\(1 + Bz — A?s*) |)? + 


const. 22) 


The constant in Eq. (22) is obtained by assuming that 
Eq. (22) must reduce to von Karman’s logarithmic 
velocity distribution law for 7/7, = 1, WV. = 
C = 0. Performing the required integration in Eq. 
(22) for these conditions and comparing the resulting 


O and 


expression with the logarithmic law we find that the 

constant for Eq. (22) must be 

const. = —1/K log, } [(T,/T.)(c,/2)]” X 
[(Pwl's)/ bel) — D (23 

Using Eq. (23) in Eq. (22) and rearranging the resulting 

expression gives 


log, } ((ewl o¥) Bol (Ty T..) (Cy 2) ] 7 = 1/1, — Dk 
(24 


ia Edr E 
o [(1 + e7)(1 + Br — A?®r?)]” (25) 


(T/T 


where 


and 
(26) 


Y(ce/2)J°E=K 


Eq. (24) can be rearranged further to obtain a compres- 
sible turbulent boundary layer velocity distribution 
with the effects of mass transfer at the wall included. 


r. 4 y\ _ 
T.. 2 6) 


Since Eq. (27) contains the local skin friction coefficient, 


This is 
. Pw. " 
6 exp (J, — DA) (27) 
Ky 


cy, and the empirical constants A and D it is of little 
The 
momentum integral including a term accounting for the 


value to us until c,, A, and D are determined. 


mass transfer at the wall can be used to determine the 
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relation between c, JJ, and R,, T,/T. and w + 
aul The constants D and A could then be deter- 
mined for the velocity distribution law only by using the 
incompressible flow Different 
values of the constant D and A will result when Eq. 


27) is used in the momentum integral to obtain skin 


velocity distribution. 


friction coefficient laws as will be seen. 


c) The Skin Friction Coefficients 

The momentum integral for the case of mass transfer 
at the wall is obtained by integrating Eq. (1) across 
the boundary layer from y = 0 to y = 6 making use of 


Eq. (7). This yields 
Cot Cy = 2(d0 dx) (28) 
where 
“l pa it y 
4 = sf ; i —- — d (5) (29) 
S Bat 2 ‘8 6 


By differentiating Eq. (27) with respect to z and sub- 
stituting the resulting expression for d(y/6) along with 
the relation between p, 7, and # given by Eggs. (11) 
and (18) into Eq. (29) we obtain an equation for cy. 
The use of this equation in the momentum integral 
implies the assumption that the turbulent velocity dis- 
tribution holds up to the wall and out to the edge of the 
boundary layer and that any error can be accounted for 
by adjusting the constants D and A such that the result- 
The 


ing expression for c, fits the experimental data. 
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cmp (~DK) 
(Cc, + c,)dx = 2d k a age | (30 
pwU WK 


MASS 


where 


me s, (2 — 2") exp (/)dz 
I, = EF : — (31 
o (1 + Cs)“ (1 + Bz —A?z2?) 


Both sides of Eq. (30) are integrated recognizing that 
the right-hand side is 0 at x = O, since E > Oas x > 0, 
and has the value of the integrand at the upper limit. 


. 


On the left-hand side C;x -f c,dx and C,x -f Cdx. 
0 0 


If 7), is not a function of x we obtain, after some re 


arrangement, 


log, [(C, + Ck. 3 + log, (K/2) — w log, (7,,/T.) + 
DK = log, (Iz) (32) 


where we assumed 

fo = told o/7 a} (33) 
Eq. (32) forms the basis for local skin friction coefficient 
and average skin friction coefficient laws. The con- 
stants D and A are determined for both the local 
and average skin friction coefficient laws by assuming 
that the appropriate form of Eq. (32) reduces to von 
Karman’'s incompressible local law* * and Schoenherr’s 


average law,* * respectively, when 7,7... = 1, M_ = 
Q and -C = 0. These laws are 
4.15 logy (c,R,) + 1.70 = cy (34) 








equation for c; is and 
Meo = 5.0 Mo = 5.0 
aa - 18 =. 1.8 
Re = 107 Re = 107 
1.0 —— T 


P coco 


| 3.19 x 10°4 
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The effect of mass injection on the turbulent boundary- 
layer velocity profile 
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The effect of mass injection on the turbulent boundary 
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Fic. 4. The effect of mass injection on the turbulent boundary- 


layer temperature distribution. 


login (C;R,) = 0.242 C, (35) 
If, in Eq. (32), we assume C, = c, and let 7,7. = 
1, \J.. = 0, and C = O and determine the constants 


K and D by comparing the resulting expression with 
Eq. (34), we find that 


K = 0.393 (36) 


and 

D = 6.53 (37) 
if we neglect terms of order 1 with respect to exp /», 
and 2/ with respect to 1, where &) = FE at 1/.. = O 
and 7,,/T.. = 1. Using Eqs. (36) and (37) in Eq. 
(32) for the /oca/ skin friction coefficient we obtain, after 
rearrangement, 


(cy + c)R. = 0.389(T,,/T .)°le (38) 


We can perform an analogous process to determine the 
constants for an average skin friction law using Eq. 
(32) with Eq. (35). In this case we find that 


K = 0.393 (39) 


and 
D = 4.13 (40) 
if we neglect terms of the same order as neglected in ob- 
taining Eqs. (36) and (37). Using Eqs. (39) and (40) 
with Eq. (32) we obtain an average skin friction law 
after some rearrangement, 
(C, + C) Re = (Tu/T.)*Is (41) 


Here J; is identical with /. except that c,in Eq. (26) is 
replaced by C, and C is equal to C,/C;. Both Eqs. 
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(38) and (41) contain the integral J, or 7; which ap- 
parently cannot be integrated in closed form rhe 
value of c, or C; for any R,, W., T,/T. and w p_U 
can be obtained from Eqs. (38) and (41) by iteration 
involving the numerical integration of J» or J. Eqs. 
(38) and (41) with /, or /; expanded into a series become 
Van Driest’s Eqs. (66) and (71), respectively, when 


mass transfer is zero.* 


(d) Heat Transfer Coefficients 


It has been shown that the assumption of an effective 
turbulent Prandtl] Number of | leads to the relation be- 
tween mean velocity and mean temperature distribu- 
tion given by Eq. (11). Using this knowledge it can be 
shown that the simple relation between skin friction and 
heat transfer known as Reynolds analogy holds in the 
present case. 

If, by definition 


du = h(T; — T,) = k,(dT/dy), (42 


where, for Prandtl Number 1, 


7, = T. {1+ [(y — 1)/2].7} (43) 
Then from Eqs. (11), (42), and (43) 
fo = h/C 2.0. =e (44) 
and, by definition 
Tw = My(du dy) x = (1/2)9..U "Cy (45) 


It can be demonstrated in a similar manner, that 
Cu = Cy 2 (46) 


Once local or average skin friction is determined using 
Eqs. (38) and (41), respectively, Eqs. (44) and (46) 
can be used to determine the appropriate heat transfer 
coefficients. 

* Van Driest’s series expansion, his Eq. (58),' must retain 
higher order terms for the correct evaluation of the integral for 
conditions of high Mach Number and large values of 7),/T 
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Fic. 6. Ratio of local skin-friction coefficient with 
transfer to local skin friction with no mass transfer vs 
transfer parameter compared with low-speed experiments.’ 
e) The Boundary Layer Thickness 


The boundary layer thickness, 6, is obtained from Eq. 


(27). Using the values of the constants previously de- 
termined for the local skin friction equation, 
0.393 and D = 6.53, the expression for 6 


namely AK = 


1S 


(6/x)R, = O.194E(T,/T..)°*! exp (h) (47) 
where here /, is the definite integral 
: ‘ dz 
I, = £ . rae (4S) 
0 (1 + Cz)(1 + Bs — A?s?)] 


The boundary layer displacement thickness, 6*, is 


defined as 
(49) 


6* is obtained using Eq. (27) in the manner similar to 
The final 
expression for the displacement thickness, using A = 
0.393 and D = 6.53, is 


that used to find the momentum thickness. 


(6*/x) Re = 0.194 (T/T 2)°*" Is (50) 
where 
i, = Ee f J ; exp (/;) . 
0 (1 + Cs)(1 + Bz — A?2s?)] 2 
b4 ) A 
cha > dz (51 
ee ‘ig )(1 oh Bz = A?2s") ( (91) 


and /, is obtained from Eq. (25). 
The boundary layer momentum thickness, 6, is ob- 
tained from Eqs. (27) and (29). For AK = 0.393 and 


D = 6.53 the expression is 
(0/x)R, = 0.194(7,/T..)*ls 


with /, obtained from Eqs. (25), (26), and (31). 


. 
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(5) ILLUSTRATION OF TURBULENT BOUNDARY LAYER 


CALCULATIONS 


Typical turbulent boundary layer mean velocity and 
mean temperature profiles have been calculated to il 
lustrate the effect of mass injection upon them. Figs 
2, 3, and 4 show velocity and temperature profiles calcu 


L.S, R, 107 


The velocity profiles 


lated for w = 0.76, M.. = 5, T,/!1 
and two values of mass injection. 


in Fig. 2 were calculated using Eqs. (27) and (3S) and 
the temperature profiles using Eqs. (27), (38), and (11 
while making the assumption that y = 6 when a = / 


in Eq. (27) to determine the constant of integration. 
While the use of Eq. (27) in this manner neglects the 
laminar sublayer, it will represent the behavior of the 
profiles as mass is injected except in the immediate 
vicinity of the wall and will serve to illustrate the effect 
It should 
that 
)]/|d(y, 6)]{ so that even if 


of the transverse mass flow on these profiles. 
be remembered in examining Figs. 2, 3, and 4 
dijdy = (U./6)} [d(a/U 
the curves shown coincided some reduction in di dy 
and dT dy would be exhibited since the boundary layer 
thickness, 6, increases with mass injection. Therefore 
a reduction in the nondimensional derivation of the 
curves shown in Figs. 2, 3, and 4 is amplified by thicken 
ing of the boundary layer into a greater reduction in 
the derivatives di/dy and dT dy at any fixed y. Since 
heat flow toward the wall and shear are proportional to 
dT/dy and di/dy, respectively, it is apparent from in- 
spection of Figs. 2, 3, and 4 that the transverse mass 
flow reduces these quantities and hence the heat trans- 
fer toward the wall and the shear stress. 

Fig. 5 illustrates the amount of mass injection neces- 
sary to reduce the local skin friction coefficient [or heat 
transfer coefficient by Eq. (44)] given amounts for w = 





2.8 } 
Moo 0 
} wos, 
: loo 
4} < 2 


9x14 < Kn < 3.3 x 10° 


SUCTION, ure a 








Cy 
7 
©Hure0 
8 | 
4} / 
| 
0 4 a 1.2 1.6 2.0 2.4 2.8 
—_ U. ea 
oo co hu 
fe Hur =O 
Fic. 7. Ratio of local heat-transfer coefficient with mass 


transfer to local heat-transfer coefficient with no mass transfer vs 
mass transfer parameter compared with low-speed experiments.’ 
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= 10’, M.,, = 5 and different values of 7,,/7 
required to 


0.76, R, 
It will be noted that the value of w/p.l’ 
reduce c, by 40 per cent is of the same order of magni- 
tude as ¢,. 

Low speed turbulent boundary layer heat transfer 
and skin friction data with and without mass injection 
have been reported in reference 9. In these tests both 
suction and injection were used at the plane boundary of 
the stream. Skin friction coefficients were determined 
from plots of the momentum thickness vs. longitudinal 
station, where the momentum thickness was deter- 
mined by obtaining the velocity profiles at various sta- 
The heat transfer was determined by direct 
Fig. 6 presents the ratio of skin fric- 


tions. 
measurement. 
tion coefficient with mass transfer to the skin friction 
coefficient without mass transfer at the same Reynolds 
Number and wall to free stream temperature ratio for 
both suction and gas injection. Fig. 7 presents analo- 
gous plots of the heat transfer coefficient ratios. The 
Reynolds Number range of these data points was about 
9 X 10% to 3.3 & 108. In all cases 7, — 7... was less 
than about 30°F. so that 7), 7, was close to unity. 
The solid lines appearing in Figs. 6 and 7 were ob- 
tained by using Eqs. (34), (38), and (44). For the case 
of W/.. = Oand 7,,/T.. > 1, Eq. (38) can be manipu- 
lated to obtain 
login [(cy + Cg) “(CyT nT) *Re| + 0.41 = 
0.242¢/((7,./T }° + (w + 1) logy 
(T./f.) Ge 


— 1[)é, 
where 
g = sin! By — sin! ap (54) 


i — 7./fJe, + &) + Td. 


a ee ’y ‘ly ‘yy ‘ly 
,(l — (T/T a) Jey — (Tu/Tw)Cq| 
and 
‘ iS C/T ep ae a of fale, ‘i 
po = —— —_ — =m (00 
ft — (7./T Jk, — (7./T cd} 
An equation for use when 7), 7. < 1 and j/.. = 0 can 
be derived in a similar manner. 
For the case of 7, 7, = 1, Eq. (53) becomes 
login [(¢, + Cy) *(cr) °R.] + 0.41 = 
(0.484,/c,) (cy + Cy) — cf?) (57) 
Eqs. (34) and (57) can be combined to obtain 
logw | [018 + @)] °{ = [0.484/(c,) “*] X 
1(1/6)(0 + o)“ — (1/2)} (58) 


where, in Eq. (58) only, 
6 = C,/ Cy; o= ¢,/t,; ¢, = 


and c,, is determined using Eq. (34). 
Because 7), 7... was greater than, but close to, one 
in the experiments reported in reference 9, Eq. (58) was 
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used in preparing Figs. 6 and 7. If Prandtl Number is 
assumed to be one, Eq. (44) applies and Eq. (58) can be 
used for Fig. 7. Reynolds Numbers of 4 X 10! and 
3.3 X 10° were used to obtain the values of c,, used in 
Eq. (58) but the resulting relations between @ and @ 
were sufficiently close that the curves drawn in can 
represent both Numbers to the 
accuracy of the figures. 
speeds, the theory of this report predicts the variation 
of turbulent boundary layer skin friction and heat 
transfer with mass transfer to good accuracy. 


Reynolds reading 


Apparently, at least at low 


(6) CONCLUSIONS 


(A) Equations are presented giving compressible 
turbulent boundary layer mean and average skin fric- 
tion and heat transfer coefficients when mass is trans 
ferred at the wall. 

(B) The equations developed are compared with low 
speed experimental results and good agreement is 
shown. 

(C) Experiments are needed at supersonic speeds to 
provide basic turbulent boundary layer skin friction and 
heat transfer data under conditions of mass transfer 
with compressibility present. 
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Buckling Stresses of Simply 
Supported Flat Rectangular Plates Under 


Combined Longitudinal Compression, 
Transverse Compression, and Shear 


JAMES H. JOHNSON, JR.* 


North American Aviation, Inc. 


ABSTRACT 


Critical combinations of longitudinal compression, transverse 
compression, and shear for simply supported flat rectangular 
platesare presented for plate aspect ratios, 1 < a/b <4 Solution 
is based on the Raleigh-Ritz method and numerical calculations 
were made by matrix iteration as discussed in reference 2. Re- 
sults presented are believed to be accurate within 1 per cent of 
exact values. The calculated data are presented in tabular form 
ind in chart form. Also included are cross plots of A, vs. A 
for various A,/A, and a/b ratios 


readily used for design and analysis. 


The latter result curves may be 


SYMBOLS 


B = symmetrical matrix |defined by Eq. (2)] ie 
D = plate flexural rigidity |Et*/12(1 — w?)], in.Ibs. 
E = Young's modulus of elasticity, lbs. per sq.in 
K = longitudinal compressive-stress coefficient 
K = transverse compressive-stress coefficient 
K = shear stress coefficient 
\ = longitudinal distributed force, lbs. per in 
NV, = transverse distributed force, Ibs. per in 
N = distributed shear intensity, Ibs. per in 
= length of plate, in. 
) = width of plate, in 


m,n, p,q = integers 
= thickness of plate, in 


plate direction coordinates 


AE = 
3 = plate aspect ratio (a/b) 
\ = characteristic value of matrix B 
= column matrix with elements x, Xe, X;, a 
a = longitudinal compressive stress, Ibs. per sq.in. 
a = transverse compressive stress, Ibs. per sq.in 
T = uniformly distributed shear stress, lbs. per sq.in 
u = Poisson's ratio 


Subscripts 


app = applied 
CT = critical 
d = design 
req'd = required 


INTRODUCTION 


Ts SKIN and other major component parts of a mis- 
sile are subjected to combined load conditions dur- 
ing flight. Buckling of parts under the load may be the 
design criteria especially when such buckling may lead 


) 
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to (1) excessive deflections, (2) redistribution of load to 
more critical members, and (3) poor aerodynamic char 
acteristics. Cover skins of supersonic surfaces, in par 
ticular, fall in this category since buckling of the skin 
may lead to poor aerodynamic characteristics and panel 
flutter. 
has been extensively treated in the literature. 


Consequently, the problem of plate buckling 


Critical combinations of longitudinal and transverse 
compression,' longitudinal compression and _ shear, 
and transverse compression and shear,” * for rectangu- 
lar plates have been found. These solutions determine 
the boundaries of the more general problem solved in 
this paper; i.e., the problem of plate buckling by com 
bined longitudinal compression, transverse compres 


sion, and shear. 


THEORETICAL SOLUTION 


The theoretical development leading to the solution 
of the general problem was made in reference 2. The 
was used with the deflected surface 


energy method 


assumed to be 


. Manx . ny 
w= Amy SIN sin 
a b 
m ln 1 


The solution is given by Eqs. (B2) and (B3) of this 





TABLE 1 
CALCULATED VALUES OF K, FOR K, = 0 
































- ° 1 2 3 | ‘ 
a*s =*sn a*on a*os | a*os } 
il 
ey even ] odd ren odd oven odd even odd | even odd 
1.0 (9.35) (12.63) 8. (6.62) | (9.61) | 4.68 | (8.42) | ) (7.05) 
1.05 438 | 
1.110 H 1899 
1.2 (8.00) | (9.70) 5.73 
1.25 7.61 | 6.0 | 1.73 4e2 
1.4286 | .% 70 
1.5 (7.07) | (7.97) 6.26 7.01 5.32 5.87 | 4.22 bods | 
1.6667 ' 3.4 33 | 
1.7 \ 5.20 } } 
1.7233 107 | 
1.75 | 6.76 6.02 | 6.19 | 
1.8 4.98 3.56 | 3.23 737 | 
2.0 (6.59) | (6.61)| 5.93 (5.72) | (5.12) | (4.66) | 4.1748 | 3.2834 | (2.89) 0 
2.2 60, 
2.3 4-43 
2.4 | 1.1255 
2.5 6.29 | (6.06) 4.6195 324273 | 3.2160] 1.22 | 1.35 
2.6 -% 
| 
2.75 | 55291 1.61 
8 4.3233 | | 
3.0 6.0% | (5.69)| 5.2775] 5.2321] 4.2420] 4.3212] 3.0167 | 3.3402 
3.2 5.0830 | | | 38579 | 
} 
304 | | 7324 . 
3.6 | | 3.09 
3.7 | | | 4ac 
4.0 5.67) | (5.77) [4.9125 |(5.00) |(<.20) |tc.22) [3.04607 | 2.9229 (2.37) 
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reference. Critical values of shear stress coefficients A, 
were calculated by allowing A, and A, to vanish singly 
from Eq. (B2) to obtain Eqs. (B5) and (BS), respec- 


Amn {(m? + n2B?)? — K,m?p? — K,n?84*] + 


where 8 = a/b; m = 1, 2, 3, aon = 1S, &......., at nt = 9D, 


In matrix form, this set of equations becomes 


where 


B 


II 


(m? — p?)(n®? — q?)V (m? + n*B? 


This system of Eqs. (1) may be solved for the lowest 
value of A, for given values of 8, A,, and A, by the 
matrix iteration method outlined in reference 2. The 
three stress coefficients thus determined define a critical 
combination of the three loadings acting simultaneously 
on the plate. 

It should be mentioned that the system of Eqs. (1) 
separates into two independent sets of equations—one 
set for m + mn equal to even numbers and one set for 
m + mn equal to odd numbers. The buckle patterns 
represented by the two sets are symmetric for m + 
even and antisymmetric for m + n odd. The lower 
value for A, from the two sets gives the critical shear 
stress in the presence of the other two loads and, at the 
same time, defines the nature of the buckle pattern. 


NUMERICAL CALCULATION AND ACCURACY 


Solution of Eqs. (1) or (2) for a given case results in an 
exact solution to the problem under consideration if an 
infinite number of equations from the system are used. 
Realistically, a finite number of equations must be con- 
sidered to obtain a numerical solution with a subsequent 


sacrifice of accuracy. Since the accuracy obtained de- 



































TABLE 2 ° 1 2 3 
—— 
CALCULATED VALUES OF K, FOR K = 0.5 a‘én a*n s*s aés 
6 
! even odd even odd even odd even odd 
= ° 1 2 3 3.5 1 8.44 7.02 5.23 2.35 
1.2 6.90 5.65 4.06 929 
| m*on a*o a*o m*on p*a 
y we odd even odd even odd even odd even odd 1.5 5.74 4.68 3.33 281 
1.0 6.73 72399 5.73 3.32 1) 1.7 ~822 
1.2 7429 6.12 4.04 2.55 
14 2.45 
1.414 0 2.0 4.97 5.73 4.14 3.11 3.54 1.43 1.59 
1.5 6.23 5.27 4.13 2.51 0,686 
1.6 1.33 2.3 3.13 3.09 
1.7 2.63 3.1 1.67 1.73 
1.8 4.6 1.95 1.27 2.5 4.73 4.98 4.00 4.07 3.14 2.91 50 
1.9 805 
2.0 5.9 6.03 | 4.8363 | 5.0597 | 4.01 36% 2.28 ° 
3.0 4.63 4.60 3.9109 3.7364 3.0131 2.6565 +22 
2.5 5.3745 | 5.36 | 4.6827 | 4.4768 | 3.87 3.45 1.89 
2.828 te) 
2.9 585 4% 3.5 4.5290 4.3872 3.7700 3.5997 2.7648 2.6084 741 762 
3.0 (5.22) |(5.04) | 4.4612 | 4.2551 | 3.5083 | 3.3228 | 2,0661 | 2.03 ° 
4.0 (4.93) | (4.89) | 4.1028 | 4.1581 | 3.1589 | 3.2936 | 1.7785 | 1.9771 ° 4.0 44132 | 4.2951 3.6244 | 3.5512 | 2.5886 | 2.6303 221 a 
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tively. If both A, and A, are retained in the de- 
velopment, the more general solution is given by the 


set of homogeneous linear equations represented by 


mnpq 


32K, 3° } } y \ 
a, = |)? 
7 . ad (> — p*)(n? — q 


n + q) are odd numbers. 


328°A,) » 


mnpq 


+ g°B")* — A,zp*B 


f 


— K,q*s* 


pends upon the finite choice made, it becomes impera- 
tive that the most important equations be chosen. A 
discussion of the guide to choosing the most important 
equations adopted by the writer is given. 

The results presented were calculated by the matrix 
iteration method using sets of sixteen equations and or 
sets of ten equations. The iterative process was carried 
through on the sixteen equations by means of the 
I.B.M. Card Punch Calculator. 
tinued until six significant figure accuracy in the charac 


The iteration was con 
teristic value of the matrix was obtained. Error in the 
corresponding value of A, was due, therefore, primarily 
to the number of equations chosen, and not to the itera 
tion. The matrix iteration method converges rapidly 
for small values of 8, A,, and A,. For these, iteration 
was made with sets of ten equations on a desk calcula 
tor. It was observed that the ten most important equa 
tions could be chosen by dropping from the sixteen 
equations those equations contributing the smallest 


T 


* Equivalent to Eqs. (B5) and (B8)?, for A; and A, equal to 


zero, respectively 





TABLE 3 


CALCULATED VALUES OF K, FOR y = 0.75 
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CRITICAL BUCKLING 


elements to the matrix. In general, this was ac- 
complished by using those combinations of m and n 
giving smallest values of the term V (m* + n’*B*) 
— K,m*B? — K,n*3* and, consequently, the correspond- 
ing term in p and g that appears in matrix B, Eq. (2). 
Values of A, obtained by the I.B.M. calculator and by 
the desk calculator were in agreement within a few 
tenths of one per cent, substantiating the method 





TABLE 4 


CALCULATED VALUES OF K, FOR y *1 


m *n = even numbers 















































| 2 ry 1 2 2.1 2.2 2.3 2.4 2.5 3 
| 1.0 8.11 | 6.63 4.685 4.19 3.63 3.33 t) 
1.1 2.45 
1.2 5.14 3.30 2.79 2.16 1.7% 
1.25 6.20 
1.3 3 
1.414 0 
1.5 5.20 | 4.00 2.25 1.66 1,28 71 
| 1,562 0 
; a .675 
| #75 4.60 
} 1.8 1.00 
| 1.826 0 
| 2.0 4.18 | 3.1476 bb 
2.1 465 
| 2.236 t) 
| 2.5 
2.6 -57% 
3.0 3.50 2.5541 839 267 
| 3.262 
} 
| 3.5 
4 3.23 | 2.3405 589%, 
oo 0 
TABLE 5 


CALCULATED VALUES OF K, FOR K = 1.5 


m #*n = even numbers 


























a, 0 1 2 
B 
1 7.39 5.73 3232 
lok 2.02 
1.1892 0 
1.2 5,56 3.91 
tes be! 
1.5 3.84 1.90 
1.6 1.07 
1.6529 0 
ig 2252 
2.0 1.45 
2.1094 0 
TABLE 6 


CALCULATED VALUES OF K, FOR Ky = 2 





1 (6.62) 4.68 
1.1 3.37 
1.2 hehh 2.01 
1.2720 0 
1.4 2.63 

1.5538 0 
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adopted for choosing the ten best equations. Further- 
more, because the use of sixteen equations increased the 
accuracy obtained only a few tenths of one per cent 
over that of ten equations, it is expected that the re- 
sults presented are accurate within one per cent of 


exact values. 


RESULTS 


Values of critical A, for all values of A, and A, for 
2S eS 
above for the symmetric and antisymmetric buckling 


$f were obtained by the method discussed 
cases. The results are given in Tables 1 through 6. 
Values shown in parentheses in the tables are values 
taken from reference 2. Values recorded to three sig- 
nificant figures were determined by desk calculator 
whereas those values recorded to five significant figures 
were determined by the I.B.M. Card Punch Calculator. 
The results presented in the tables are repeated in 
graph form in Fig. 1. The conditions for symmetric 
and antisymmetric buckle patterns are clearly shown. 
For brevity, regions of symmetric buckle patterns are 
designated by 1J = 1 and J/ = 3. 
symmetric buckle patterns are designated by .1/ 2 


Regions of anti- 





FIG. 1(a) CRITICAL STRESS COEFFICIENTS FOR FLAT 
RECTANGULAR PLATES UNDER COMBINED LOADS 
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FIG.1(b) CRITICAL STRESS COEFFICIENTS FOR FLAT 
RECTANGULAR PLATES UNDER COMBINED LOADS 
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FIG. I(c) CRITICAL STRESS COEFFICIENTS FOR FLAT 
RECTANGULAR PLATES UNDER COMBINED LOADS 
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FIG. 1(d) CRITICAL STRESS COEFFICIENTS FOR FLAT 
RECTANGULAR PLATES UNDER COMBINED LOADS 
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FIG. I(e) CRITICAL STRESS COEFFICIENTS FOR FLAT 
RECTANGULAR PLATES UNDER COMBINED LOADS 
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— (5) Determine margin of safety. Assuming that th 
FIG. 2(¢) DESIGN CURVES FOR THE BUCKLING OF FLAT wi 
RECTANGULAR PLATES UNDER COMBINED LOADS loads increase at the same rate and are therefore in th« 
o/b +10 same proportion to each other at all load levels, the 
J margin of safety based on load is given by 
M.S. = (Ne/Napp.) — 1 = (ta /treara l 
Margin of safety based on stress is given by 
al M.S. = (Gce,/Capp.) — 1 = (ta/trega)? — 1 
cr./ Tapp d/ req'd 
A NUMERICAL EXAMPLE 
» bo 
N_ = 32 Ib/in 
K, Y 
| 
and .1J = 4. Furthermore, the curves, as well as the ! 
tables, show that for A, > 1.0 only symmetric buckling —e} | N = 80 |b/n om 
is possible for any value of a/b. These curves may be N. = 100 Iba) | ” | ss 
used to determine critical combinations of stresses not | * ol 5 
given by the table. —_— | | 
Although the curves of Fig. 1 show the trend of sym- ~ _ ™ f 
metric and antisymmetric buckling to good advantage 
and permit interpolation of results between those | 
determined by calculation, they do not present the data re a a 


in a manner convenient for use in design. A plate in a 
critical state (neutral stability) under combined loads 
must carry stresses in direct ratio to the applied loads; 
tet as, NN, = B/K, and B/N, = &,/K, A 
critical combination of loads is defined by the point:on 
the curves corresponding to the ratios A,, A, and K,/K, 
at a given value of a/b. 

A convenient way of plotting the data to easily deter- 
mine the critical combination is used in Fig. 2._ Data for 
these figures were obtained by cross-plotting the data of 
the tables and Fig. 1. Sufficient points have been cal- 
culated to allow the cross plots. To use the curves, the 
following procedure is to be followed: 

(1) Calculate the ratios K,/K, = N,/N, 
K, = N,/N.. 

(2) On the curve corresponding to the plate a 0, lay 
off a straight line from the origin of slope A, A,. 


and A, + 


(3) At the intersection of this line and the curve 
corresponding to the A,,/A, ratio of Step (1), read A, 
and/or A,. 


(4) Determine required plate thickness from 


’ rE i 
N ; = (Or) app.t = Zz ° 9 
12(1 — p*)b 
or 
} = (Nz)app 12( l — u*)b° ae (Ny) app 12( C= u”)b? = 
_— K, rE K, rE 


(N.)app.12(1 — u2)b? 


K, rE 


The value of A, may be determined from A, and K,/K 


if desired. 


Consider the plate shown with four simply supported 
edges anda = 10in., 6 = 5in., ¢ = 0.051 in., uw = 0.30, 
E = 10° Ibs. per sq. in. Determine the plate margin 
of safety. 

Calculate the stress ratios and load ratios 

K,/K, = N,/N,; = 32/80 = 0.4 
K,/K, = N,/N; = 100/80 = 1.25 


On Fig. 2c, the interaction curves for a b = 2, lay off 
a line from the origin of slope A,/A, = 1.25. The in 
tersection of this line with the curve A, A, = 0.4 de- 


termines the critical buckling coefficients for the three 


loads. From Fig. 2c, 


R= 25: A, = 20: and K,/K, = 0.4 


therefore K, = 0.8. 





FIG 2(b) DESIGN CURVES FOR THE BUCKLING OF FLAT 
RECTANGULAR PLATES UNDER COMBINED LOADS 
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To determine the plate thickness required to sustain 
these loads, any one of the three buckling coefficients 
determined may be used. Using K,, 

. rt 
t{=K.- —— 
12(1 — p?)d? 


or 
N,12(1 — p?)b?]'? : 
are 7 7 | = 0.048 in. 
ie rk 


Since the actual design thickness is 0.051 in., the margin 
of safety based on stress is 


The results presented herein are believed to be ac- 
curate within one per cent of exact values. 

Interaction curves of Figs. 2 may be used to deter- 
mine critical combinations of shear and direct stress for 
buckling of flat rectangular plates. It is recommended 
that the method outlined for the use of the curves be 
followed for design and analysis. 
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Differentiation of Experimental Data 


Peter F. Jordan 
The Glenn L. Martin Company, Baltimore 3, Md. 
February 1, 1954 


mgt of experimental data is the subject of re- 
cent contributions by Kauffman and Shinbrot! and by 
Wand and DeSanto.? 
The second authors give equations in terms 


The first authors propose to use a semi 
graphical method 
of higher differences; they give equations both for inner points 
and for end points of the given experimental curve. The present 
remark concerns only the case of inner points. 

The writer tends to agree with the second authors in that 
graphical interpolation would seem to be a little out of fashion in 
the present age of electrical desk machines; on the other hand, 
he contends that higher differences, though they are an excellent 
means for writing down general equations, are not liked either by 
What he likes to use are 

Reducing both proposed methods 


the practical calculator sums of 
weighted functional values 
to this form, the first by means of mechanical interpolation, 
the second by dissolving the higher differences, we find that the 
two methods are in fact equivalent in that both lead to the same 


result, viz., to 


f"(X6) f’ (Xo (1/12h)(f_2 — 8f.1+0+ 8fi —f (la) 
"un ag € (1/12h?)( —f_» + 16f_, — 30f) + 

16f, — f (1b 
Here 


- f(x» + nh 


As had to be expected in this well-covered field, the relations (1 
are by no means new. The complete equations, see, e.g., 


Schulz,* p. 128, read: 


f'(xo) = f’(xo) + (h*/30)m, 4 / 
(2) 
f"(x0) = f"(x0) + (h*/90)me 4 \ 
where 
mr| S (max.|f''7(X)| xg —2h < x < xo + 2h 


From Eq. (2) it appears that Eq. (la), using only four func 
tional values, is still correct for 4th order polynomials. The 
explanation for this advantageous property of Eq. (la) which, 
in combination with its simplicity, makes Eq. (la) well suited for 
practical use, is that f’(x,) = 0 for any 4th order polynomial with 
= f-, =f, = fe. A corresponding argument applies to Eq 

lb) which formula is still correct for 5th order polynomials. 
Considering the task of differentiating a given experimental 
curve by means of Eqs. (1) it will be noted that the effect of a 
given inaccuracy A of one of the functional values f, decreases if 
his increased; at the same time the error term in Eq. (2) in- 


creases. The optimum value for A can be estimated if both the 


order of A and the nature of the function f(x) are known 
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The Flow on a Rotating Body of Revolution in 
Axial Flow 


H. Schlichting and E. Truckenbrodt 

Institute of Fluid Mechanics, Technical University, 
Germany 

January 22, 1954 


Braunschweig 


6 ip FLOW IN THE BOUNDARY layer of a rotating body of revo 
lution including the rotating disc in an axial main flow has 
been calculated by the authors both for laminar and turbulent 
incompressible flow in several German publications. A summary 
of the results is given here 

For the rotating disc the solution is valid for 0 Veal Ua * : 
whereas for the body of revolution the solution is confined to 
values of V,,/U., not too large ( V,, circumferential velocity of 
the main cross section, U’,, axial velocity), see Fig. 1 

The rotating disc in axial flow and with a laminar boundary 
layer has been dealt with by an approximate method by H 
Schlichting and E. Truckenbrodt 
problem was published by D. M 
and Sheng To Chu 
The extension of the solution for the disc to turbulent 


Later an exact solution of this 
Hannah? and A. N. Tifford 
There is a good agreement of these two 
solutions 
flow is due to E. Truckenbrodt 
is basically the Pohlhausen momentum theory, which in the 


The approximate method used 


turbulent case requires some further assumptions about the 
shearing stress at the wall 

The main results are explicit formulas for the moment coeffi 
cient Cy = M/(p/2)w?R® (w = angular velocity of the disc, R = 


radius of the disc One obtains for the disc, wetted on one side 


only, (V»/U. < 1): 

Cy = (2.53/V Re,)(U./Vm)*5 laminar (1 
and 
(turbulent (2 


Cu = (0.079/W Re) Us, / Vin) 


Here V,, = Rw and Re, = R*w/v is the Reynolds Number 
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Fic. 1. Notations. (a) disc; (b) body of revolution. 

The rotating body of revolution with the laminar boundary 
layer has been treated on similar lines by H. Schlichting.6 One 
obtains two momentum equations for the meridional and circum- 
ferential direction. By introducing polynomials for the two 
velocity profiles with appropriate boundary conditions one finally 
arrives at two simultaneous ordinary differential equations which 
are solved numerically for different values of V,,/U,,.. The main 
results are the skin drag coefficient and the moment coefficient, 
as they depend on V,,/U.,, and the Reynolds Number. As an 
example the rotating sphere in a parallel flow along the axis 
of rotation has been worked out, and the moment coefficient is 


found 
M 815 Uv. 


Cu = = (laminar ) 
V Re V» 


(p/2)w*R® 


(3) 


with Re = U,.R/v; Re, = R*w/v; Vm = Rwand R = radius of 
the sphere. * 

The extension to the rotating body of revolution in axial flow 
with a turbulent boundary layer has been given by E. Trucken- 
brodt.?. The turbulent case again requires some additional as- 
sumptions about the shearing stress at the wall for both the merid- 
ional and circumferential component. These are deduced from 
an analogy with the yawing wing of infinite span. 
E. Truckenbrodt succeeded in simplifying the calculation to such 


an extent, that the integration of the two simultaneous differen- 


Furthermore, 


tial equations mentioned above is reduced to the following quad- 


rature formula for the moment coefficient Cy = 17/(p/2)V,,?R,,2/ 
for an arbitrary body of revolution (R,, = radius of maximum 
cross section, / = length of reference, V,, = R,w). 


* The rotating sphere in a fluid at rest has been dealt with recently by 


I. Howarth. Our method is not applicable to this case 
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Here R(x) is given by geometry, U(x) is the potential flow for the 
is the skin friction coefficient of the flat 
plate for turbulent flow at Reynolds Number U.,//v. Th 
constant KA, stands for the contribution of the laminar part 


body of revolution, Cc; 


O<x < x, given by 


; jl “aul oy (R\6 x ]/2y) 
Ky = < Cri : ) d 
\2 R l R, i} f 


where cy; is the laminar skin friction coefficient at Reynold 
Number U,//v. For laminar flow the exponent is equal to 1, 
whereas for turbulent flow » = 1/6 
has to be taken until the trailing edge x» 

Eq. (4) is valid for laminar and for turbulent flow as well, and 


In Eq. (4) the quadratur 


is similar to a quadrature formula for the momentum thickness 
for a nonrotating body of revolution, due to E. Truckenbrodt. 
It contains also, for x — 0, the rotating disc in axial flow for mod- 
erate values of V,,/l For the rotating sphere Eq. (4) is in 
agreement with Eq. (3) for laminar flow, whereas for turbulent 
flow one obtains with / = R and Re = U,.,,R/v and with the 


same definition for the moment coefficient as in Eq. (8 


7 , , 
Cu = (0.299/V Re)(L Vim) (turbulent 6 
As further examples rotating ellipsoids have been worked out, 


with the major axis as axis of rotation. In Fig. 2 the moment 


coefficient Cy = .M/(p/2)V,?R,?L is plotted against Reynolds 
Number Re = U.L/v for different values of k = D/L. Here 
V, is the circumferential velocity of the minor axis, V,, = R,,w 


and L the length of the body. 


REFERENCES 
1 Schlichting, H. and Truckenbrodt, E., Die Strémung an einer angeblasene 
votierenden Scheibe, Zeitschrift angewandte Mathematik und Mechanik 
Vol. 32, p. 97, 1952; Journal of the Aeronautical Sciences, Vol. 18, p. 638 
1951 
2 Hannah, D. M., Forced low against a rotating disc, A.R.C. Rep. a. Mem 
2772, 1952. 
Tifford, A. N. and Chu, Sheng To, On the Flow Around a Rotating Disc i 
a Uniform Stream, Journal of the Aeronautical Sciences, Vol. 19, p. 284, 1952 
* Truckenbrodt, E., Die Strémung an einer angeblasenen rotierenden Scheibe 
bei Lurbulenter Strémung. To be published in Zeitschrift angew. Math. Mech 
1953 
5 Schlichting, H., Die laminare Strémung 
tierenden Drehkor pern, Ingenieur-Archiv, Vol. 21, p. 227, 1953 
given at the VIII. International Kongress of Applied Mechanics, Istanbul 


um einen axial angestrémte 


From a paper 


1952 
© Howarth, L., 
VII, 42, p. 1303, 1951 
Truckenbrodt, E., Ein Quadraturverfahren zur Berechnung der Reibung 


Note on the boundary layer on a rotating sphere, Phil. Mag 


schicht an axial angestrémten rotierenden, Drehkorpern, Ingenieur-Archiv, 
vol. 22, p. 21, 1953 

8 Truckenbrodt, E., Ein Quadraturverfahren cur Berechnung der | 
Reibungsschicht bei 


Ingenieur-Archiv, Vol. 20, p. 211, 1952; 


amunarent 
eben:r und rotationssymmetrische 


Journal of the Aero 


und turbulenten 
Sirémung, 
nautical Sciences, Vol. 19, p. 428, 1952 


The Compressible Laminar Boundary Layer 


Deane N. Morris 

Aerodynamics Engineer, Douglas Aircraft Company, El Segundo 
Calif. 

January 28, 1954 


IT OUR RECENT PAPER on the compressible laminar boundary 
layer,! the numerical results had to be restricted to the use ofa 
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READERS’ 


linear Viscosity law, although the theoretical analysis did not re 
quire this restriction. By employing a stagnation temperature 
profile instead of a static temperature profile in the analysis, the 


iculties with regard to the use of an arbitrary viscosity law 


diff 
which were discussed in the reference paper) have now been 
eliminated. Apparently the stagnation temperature, being 
much more slowly varying through the boundary layer, is more 
idaptable to the type of analysis of reference 1 than is the static 
temperature 

In line with current work on the laminar boundary layer, the 
theory also has been extended to allow arbitrary variation of 
specific heat, thermal conductivity, and Prandtl] Number, and 
the use of imperfect gas relations 
omplications in the work, requiring only the use of stagnation 


These changes cause no basic 


enthalpy in place of stagnation temperature and certain addi- 
tional terms of the type employed for the arbitrary viscosity law 

A report is in preparation which uses stagnation conditions for 
the basic thermal variable and includes the variable properties 
listed above 

In reference 1, values of g,(0) and g(0) had been obtained as 
igenvalues of the equations representing the initial conditions 
for finite leading-edge velocity and a stagnation point, respec 


tively. It is now clear that the equation for the stagnation point 


does not possess a singularity as 7,,* — 1 and um — 0 (as stated in 


reference 1) since mu at x = O is always identically zero. There 
fore no special conditions need be imposed upon g.(0), and it 
issumes the constant value of —1/60 for a sixth-order poly 


nomial 
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Some Interference Effects Between Two Flat 
Surfaces Planing Parallel to Each Other at 
High Speed* 


Daniel Savitsky and David A. Dingee 

Experimental Towing Tank, Stevens Institute of Technology, Hoboken, 
N.J 

January 15, 1954 


- THE COURSE of recent experimental planing investigations 
conducted at the Experimental Towing Tank, some data 
were collected on the change in lift and wetted area experienced 
by a flat planing surface when operating parallel to, and at various 
lateral distances from another similar surface at identical planing 
conditions. These data were collected over a limited range of 
test conditions and hence are not inten ‘ed to represent a compre- 
hensive study of the interference effects between parallel-running 
planing surfaces. Nevertheless, it is believed that a presentation 
of the results will be useful to indicate the type and order of 
magnitude of the effects to be expected, and the lateral spacing 
between surfaces for which these effects are important 

A similar interference phenomenon has been considered by the 
1erodynamicist in the case of two or more wings of high aspect 
ratio in flight side by side. It can be demonstrated that a result- 
ant upwash velocity is induced by the trailing vortices of one 
wing on the other giving each wing a larger effective incidence 
angle, relative to the free stream, than would exist if the two wings 
This effect naturally results in an 
The curve of Fig. 19 of 


were infinitely far apart. 
increased lift efficiency for each wing 


* The material reported herein was obtained in the course of the Research 
on Planing Surfaces at the Experimental Towing Tank of Stevens Institute 
of Technology under contract with the Office of Naval Research 
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TABLE | 


Test Conditions 


ce T h Ait 
7.50 8 0.325 2 60 
7. 50 8 0.596 4 23 
7.50 12 0.256 1.20 
7.50 12 0. 504 2 40 
*( V/V gb 
TA h/b sin + 


reference 1 implies that there will be a very rapid increase in lift 
as the lateral spacing between two high aspect ratio wings is 
reduced. An analogous analytical study of this phenomenon in 
the case of planing surfaces is complicated since the craft is a low 
aspect ratio lifting surface which operates on a free water surface 
whose shape is deformed by the interaction effects 


Tests were made on two 4-in.-beam flat surfaces at the planing 
conditions given in Table 1, with the lateral distances (d) between 
surfaces varying from 0 to 5 beams. In order to establish an 


end poit a single planing surface was also tested 


Since the subject interference effects are due to the dynamic 
component of the planing force, the tests were made at a relatively 
high speed coefficient (C, = 7.50), where the buoyant contribu 
tion to lift is small. The test procedure was to measure the lift 
and photograph the wetted bottom area of a single planing surface 
running alone at a fixed heave and then, keeping all conditions 
the same, to measure the lift and wetted area when another iden 
tical surface was located parallel to and at various lateral dis 


tances from the first surface 


INTERFERENCE EFFECT ON MEAN WETTED LENGTH-BEAM RATIO 


In the case of flat planing surfaces, the mean wetted length 
beam ratio (A) is defined as the distance, in beams, from the trail 
ing edge to the stagnation line intersection with the bottom, (A 
is the inverse of the aerodynamic definition of aspect ratio AR 
An examination of mean wetted length data obtained in tests by 
Sottorf, Sambraus, and in recent NACA and ETT tests has 
indicated that, for flat planing surfaces, with no interference 
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effects, the running mean wetted length-beam ratio (A) is 0.30 
larger than the mean wetted length-beam ratio measured to the 
level water intersection with the planing surface (A,). This so- 
called wave rise relation, 68 = A — A; = 0.30, has been found to 
be generally applicable for \; > 1 over a wide range of trim angles 
and speed coefficients and includes the test conditions reported 
herein. The disturbing effects of the upwash velocities, resulting 
from the interference phenomena, on the wave rise can be seen in 
Fig. 1, where the quantity 6A is plotted against separation dis- 
In this plot, the beam (}) of one surface is used as the 
The data for the test trims and 


tance. 
nondimensionalizing factor. 
wetted lengths collapse on a single curve fairly well and indicate 
a wave rise variation starting from a value 6A = 0.30 at about 
three beams spacing and increasing to 6A = 0.60 at zero beam 
spacing. At this zero spacing end point, the wave rise can be 
shown to be in agreement with the original relation 6A = 0.30 
if, in this case, the wave rise is expressed in terms of the geo- 
metric beam of the combination of surfaces, 2d. 


INTERFERENCE EFFECT ON Lift COEFFICIENT 


A plot illustrating the experimental variation in lift with 
lateral distance between two flat surfaces planing side by side is 
given in Fig. 2. The ordinate of this plot is the ratio of the actual 
measured lift (A,) for one flat plate running parallel to another 
flat plate at identical conditions of trim, speed, and depth of 
immersion, to the lift (A) for one surface planing alone at the 
same trim, speed, and depth of immersion. It is seen that A; 
is always greater than A. This fact may be attributed, essen- 
tially, to two interference phenomena—one is the influence on lift 
of the change in upwash velocity, while the second is the influence 
on lift which results from the wetted length change, as discussed 


in the previous section. 


The collected test data collapse on a single curve very well and, 
similar to the wetted length variation, there is a lift increase on 
each surface, due to proximity effects, for separation distances 
between zero and three beams. This lift increase, as the spacing 
is reduced, is not nearly so rapid as the corresponding lift increase 
for high aspect ratio wings. At spacings greater than three 
beams between surfaces, the interaction effects are negligible. 
When both surfaces are in contact with each other (zero beam 
spacing), the test data show the lift increase to be approximately 
47 per cent. This end point is in agreement with the results 
that would be obtained by substituting the present test param- 
eters into empirical planing lift equations such as those given in 
references 2 and 3. For values of \ outside the present test range 
some change in the lift interference curve of Fig. 2 is expected. 
The ratio A,/A at the end point d = 0 can be calculated for these 
other values of \ by using the aforementioned empirical equa- 


tions. 


AERONAUTICAL 
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It is not expected that the curve of Fig. 2 will be applicable for 
low speed coefficients, where the buoyant contribution to lift js 
appreciable, since, at zero velocity, where the entire load js 
supported by buoyant forces, there is no variation of lift with 
lateral spacing between the surfaces. 

As stated above, one contribution to the lift on a planing body 
in proximity to another results from a wetted length change. It 
was found that, for the conditions tested, consideration of the 
effect of increased wetted length alone could account for approxi- 
mately 15 per cent of the total lift increase at any lateral spacing 
of the surfaces. Thus, as in the case of wings, the charge ijn 
downwash velocities is of primary importance in affecting the lift 


on planing surfaces operating side by side. 
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Methods for Determining the First and Second 
Derivatives of Experimental Data 


Ernest E. H. Schurmann and Theron L. Smith 

Aerodynamics Engineer and Aerophysics Engineer, Respectively, 
Consolidated Vultee Aircraft Corporation, Fort Worth Division, 
Forth Worth, Texas 


January 22, 1954 


I TWO EARLIER NOTES in the Readers’ Forum, June, 1953, and 
November, 1953, two methods for the differentiation of ex- 
perimental data were presented. The first note, by Kauffman 
and Shinbrot,! described a method which used first order differ- 
ences and a graphical interpolation to calculate the derivatives. 
The second note, by Wang and DeSanto,? found that by using 
higher order differences some of the shortcomings of the first 
note could be overcome. 

Upon investigation of both methods it was found that if a 
third degree equation was used to fit the data at the four points 
being used, the method presented by Kauffman and Shinbrot 
reduced to that used by Wang and DeSanto. As was pointed 
out by Wang and DeSanto in order to evaluate the derivatives 
at the front and back portions of the curve, a method of forward 
and backward differences must be used. 

The authors found that the equations used by Wang and 
DeSanto reduced to those given in reference 3. The equations 
given in reference 3 are in form that makes them easier to handle. 
Wang and DeSanto use a single equation to calculate the deriv- 
atives at the first and second points while, in reference 3, two dif- 
ferent equations are used. The values calculated at the first 
point agree, but much better agreement can be obtained at the 
second point if the equation from reference 3 is used.* The same 
is true for the last and second from last points 

The equations that are given in reference 3 are reproduced here. 
The notation used is shown in Fig. 1 
(dx/dt)att, = 1/24 At (Ax_2 + Bus + Cxo + Dui + Ex2) (1) 


(d2x /dt?)att, = 
1/[12(At)*] (Ax_2 + Bry + Cxo + Dur + Ex: 


* Using method in reference 2, the authors could not check the answers 
calculated by Wang and DeSanto for the first two and last two points on the 


dx/dt and d?x/dt? curves 
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The coefficients that are used in these equations are given in 


sable for TABLE | zs tae sein 
to lift is Cuiiivates ter ta, 1 Tables 1 and 2 The above form for the derivatives can be easily 
load is adapted for a high-speed electronic digital computer 
ift with n A B eo s D E The method was checked by using the same problem that was 
= —50 96 _79 39 _6 used in references 1 and 2, namely, 
ng body —1 —6 —20 36 —12 2 x = 6sin¢ + 0.5 sin 3¢ (O<t<s3 (3 
NO 9» an > > an 
a Bs ; _9 _ uae = 6 pncrements of 0.2 int were used and the input data was used cor- 
ation 9 6 39 > _96 50 rect to 3 decimal places. rhe results are shown in Fig. 2 
4 ae As is shown in Fig. 2, the results are in excellent agreement 
alten. with the theoretical curve even at the ends of the curve. In the 
4 - middle portion of the curve the root mean square errors were 0.005 
the lift TABLE 2 and 0.037 in dx/dt and d*x/dt?, respectively. At the ends of the 
Coefficients for Eq. 2 curve where backward and forward differences were used, the 
= root mean square errors were 0.004 and 0.120 for dx/dt and 
n A B ( D E d?x /dt?, respectively 
pe —2 35 —104 114 — 56 11 
liam F.. an 11 2) 6 { = REFERENCES 
titute of 0 —] 16 — 30 16 —] 1 Kauffman, W. M., and Shinbrot, M 4 Method for Differentiation of 
st 1949; 1 a 4 6 — 2) 11 Experimental Data, Journal of the Aeronautical Sciences, Vol. 20, No. 6, pp 
oo 2 11 —56 114 —104 35 Se ee 
: 2 Wang, Chi-Teh, and DeSanto, D. F., Differentiation of Experimental Data 
Lifting by Means of Higher Order Finit.-Difference Formulas, Journal of the Aero 


nautical Sciences, Voi. 20, No. 11, pp. 792-793, November, 1953 
3 Southwell, R. V., Relaxation Methods in Theoretical Physi Ist Ed., p 


230: Oxford University Press, London, 1946 
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| | Load Distributions on Some Twisted and 
x lx Ix Ix lx Cambered Delta Wings with Supersonic 
tivel - * fe , e Leadi d Trailing Ed 
vely, | | | eading an railing ges 
| 
VISION, | | | At 
| 
1 1 , | t H. K. Zienkiewicz 
to t, to t, to Aerodynamicist, English Electric Co. Ltd., Luton, England 
” ; January 19, 1954 
», and Fic. 1 
of ex- ‘ 
ff SUMMARY 
man 
differ- 8 q OA calculated Using Puckett’s source distribution method, integral relations are de 
tive veloped for the load distribution outside the tip and subsonic leading edge 
— — true regions of symmetrically twisted and cambered wings with supersonic leading 
using ie snag edges. These relations are used to calculate the load distribution on delta 
first ‘\ wings with the incidence distribution proportional to x” y%for0 < p,q S 3 
dx O0<p+a53 
t if a \ . 
. t NOTATION 
01Ints 4 
nbrot Ap/@ load (lifting pressure) coefficient 
inted f(y), a(x incidence distribution functions 
a F(y = Pg. y) dy 
AtIVES H(é,, x,y) = [(x £)? B2(y — n) 
‘ward i Hi(n, x,v) = {la — LO)? 32(y r 2 
Ae, x, ¥ i(x — &)? — p?ly I €)}? 
h arbitrary constant 
and : z ; 
k Ptr ¢ 
tions O* je functions defined by x Liy) and y I x) along 
J 
ndle. 30 paths of integration (Fig. 1 
leriv- \ x j n cot 8) /B 
i \ 2 p,q = non-negative integers 
o dif- \ dt P, 
. \ t = By/x 
first " z = 2(x, y) equation of wing surface 
> \ } 
t the ' a 8 = VM?-1 
same 4 1 0 = wing apex semiangle 
d r§ / r 02/0x, local wing incidence 
here. i * b T ‘ ; ¢ = velocity potential on wing upper surface 
\ / 4 / ¥ = cos! [(m + t)/(1 + nt)] 
\ / \ / 
\ / \ 1 Primes denote differentiation with respect to the argument 
\ 
( l \ rat "4 4 
\ 
; 4 5 ANALYSIS 
® . 
2 \ ‘ \ / (1) General 
; ‘ 
\ : F ; ; 
i wl — HAS SHOWN! that the velocity potential, ¢, at a point 
poe -8 P(x, y) on the upper surface of a thin wing (Fig. la) whose 
n the > e - . F 
Fic. 2 leading edge is supersonic between A and B is 
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Fic. 1. Paths of integration for various incidence distributions 


y=L-\x). (b)A=hly 


¥ — 
g= J / AE, n)H(E, n, x, y) dé dn (1) 
T Ss 


where S is the region AOBP of Fig. 1. The load coefficient is then 


Ap/q = (4/V) (0¢/ox (2 
As shown by Evvard,? 0¢/0x can be expressed in the form 
mo rr V “OA 
= hHdn + H —dé dy 3) 
Ox 7 JAOB 7 Ss og 


where — and X are evaluated in the line integral as functions of 
n along the path of integration AOB. 

We shall now consider certain types of the local incidence dis- 
tribution, A, and give expressions for Ap/g which involve line 
integrals only. 

(i) X = hf(y); fly) = f(—y). 
rically twisted wing and since OA /Oé 


This is the case of a symmet- 
= 0, the surface integral in 
Eq. (3) vanishes and we have, simply, 
am Ap , 
= f(n)Hi(y, x, y) dn (4 
4h q AOB 


(ii) X = hg(x).—In this case the wing is cambered but not 


twisted. Substituting for \ in Eq. (1) integrating with respect 


to 7 and differentiating partially with respect to x, it can be shown 


that 
ra ~ fy — L-\E)] 
p = | = z g(£)H2(E, x, y) dé (5) 
4h q AOB oe 
(iii) X = hxf(y); f(y) = f(—y).—Such wings are twisted and 


cambered. Substituting for \ in Eq. (3) and integrating with 


respect to ¢in the surface integral we find that 
a Ap 


Lin) f(n)Hi(n, x, vy) dn + 
a. n)f()f1i(n y) dn 


AOB 
F(n) 


[Liy) — x] J 
AOBY~— 1 


Comparing Eqs. (6) and (4) it will be seen that Eq 


Hy(n, x, y) dn (6) 


(6) can be 
written as 
Ap A 
P=? iLyyty) + 
q q 


[L(y) — x] - (6a) 


th F(n) 

H dn 
= AOB ¥Y —7 
where (Ap/q) [L(y)f(y)] is the load coefficient for a purely twisted 
wing with the incidence distribution \ = hL(y)f(y). 





~ 47 








Poy) 
~% 


(a) X = Af(y), X = hg(x), and A = hxf(y 
g(x); along OC, L-\x) = 0 


along AOB, 


(iv) A =h 
bined twist and camber 


This is another case of wings with com 
With this A, Eq. (3 


¥) gl) 
leads to 


t 


a Ap nw Ab 
= — tly) g(Liy)]}} 


bens ( 


4h q 4h q 


7 ~ “ 
ge Qs | 

yly — L-\(x)] | -— Hodé 4 / Hdé | - 
CcCOAy “Ss con + ¢ J 


~ + 


1 = 
l J g'(é)Hodt — — g'(E)Hed= + YUN 7 
8? Jooa B° Joos 4 


where Ap/g}| yi g{L(y)]} is the load coefficient for a purel; 
twisted wing with the incidence distribution \ = h y g{L(y 
and y = L (x) is the equation of the paths of integration COA 


and COB (see Fig. 1b 
(2) Delta Wings with Supersonic Leading and Trailing Edges 

The general relations developed above have been used to find 
the load distributions on symmetrical delta wings with super 
sonic leading and trailing edges and with the incidence distribu 
tions of the form \ = hx?|y/4, forO < pg 53,0 Sk = 
g <3. The results show that for \ distributions of this type the 
load distributions due to camber and twist can be expressed 


J(t??) VQ1l — n?)(1— #@)4+ A(f#)sech'/t); O< #2 <1 s 
1<f < a? 


(2/4h)(B/x)(Ap/q) = G(t)y(t) + G(—tyYy(—t) + j 


= rG(—'t!); 


G, J, and K are algebraic functions of ¢, n, and 8, tabulated in 
Tables 1 and 2. Our result for the wing with linear camber 
(p = 1, q = O) agrees with the well-known solution obtained by 
several investigators (e. g., Heaslet and Lomax’). Since the 
work reported here had been completed, Kainer* has published 
his solutions for wings with symmetric linear, quadratic, and 
cubic twist (bp = 0, q = 1, 2, and 3). Our results agree with 
Kainer’s except that the load coefficients as given by Fq. (8 
are 1/8 times those given by Kainer. Now for the limiting cas¢ 
of flat rectangular wing (the region unaffected by the tips 
Kainer’s results reduce to Ap/q = 4a, where a@ is the wing inci- 
dence, while the correct result is, of course, Ap/g = 4a/8. Thus, 
it appears that our results are correct and that the left-hand sides 
of Eqs. (11), (13), and (17) 


Figs. 2 to 9 of reference 4 should be multiplied by 6 


and the corresponding captions of 


An interesting point arising from the present investigation is 
that, out of all the cases considered here, the sech™!| ¢) term ap 
pears in Eq. (8) only for the incidence distributions proportional 
and it can be shown that, in general, the fun 


tox|y| and x?/y 
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PABLI 


The Function G 


1 — n*) Git 
B 
-(1/8 (% + ¢ 
— 2n? — n*t 
2 1/28) [2(n + t)? + (1 + nt 
1/2) [2n(n + t)? — 2(1 — m?)?t(1 + mt) + 


n*(1 + nt 


B/2) [B8n2(1 + nt)? + 2(1 — n?)? — 4n7%(1 — mn?) & 
l + nt 
—(1/2g n+t) [2(n +t)? + 301 + at)?] 
2 1/6) } —6n(n + t)3 + (1 + ant) [601 — n?)2t(n + 4) - 
Qn(1 + nt) (n + t) + (1 — n?) (2n? + 1 1 + mt)?]j 
9 —(B/2)|) }2n%n + t)? + 4(1 — m?)3(1 + ant 
(5n2 — 2)t + n¥(1 + 2n?)] (1 + nt)?} 
32/2) 211 — n?)® — n%1 + nt) [601 — n? _ 
ai] — n? l nt) + (n* + 4 1 + nt)?]; 


TABLE 2 


The Functions J(f?) and A(?? 
l1—*» 2) I(t? K(t 
() () () 
] 2 3 () 
a an 0 
2 —(32/B () 
1 1 — 4n? t 
9 () Bn(2 — 5n? 0 
; 1/38) [4 + lln? + (4n? + 11)P 0 
2 1/3)n|] [1 + 140? + (6 + 11 n? — 2n‘)t? 0 
2 1 (8/3) [2 — 6n? + 19n1 — (2 — 9n? — 8n‘*)t? 2Bt? 
Q (B2n/3) [n*(2 + 13n?)t? + (6 — 17n? 26n‘*)) 0 
tion A vanishes except when A=x? y'*4*1,g >0,p>0. Thisis 


in agreement with our results? for symmetrically twisted and 
cambered delta wings with subsonic leading edges, where a 
similar state of affairs is found to exist; in the limiting cases of 


sonic leading edge (wn = 1) the sech™! |¢ terms obtained from 
the supersonic leading edge case become identical, as they should, 
with the corresponding terms obtained from the subsonic leading 
edge case This provides a further check on both sets of 


results 


(3) Applications to Other Wing Plan Forms 


The above results also hold, of course, outside the tip and sub 
sonic trailing edge regions of sweptback wings with supersonic 
leading edges 

The limiting case of a wing with unswept leading edge (regions 
outside the influence of the tips) can be readily obtained by setting 
n = (in Eq. (8) 
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Flow Between Nozzles in Series with an 
Isentropic Work Change 


Charles A. MacGregor 
Aerodynamic 

Calitornia 
January 19, 1954 


SYMBOLS 


i irea 

( specific heat at constant pressure 
absolute total pressure 

R gas constant 

1 absolute total temperature 

W weight flow of gas 


ratio of specific heats 


difference between two value 


“be 


ratio of per cent of sonic flow 
Subscripts 


critical (sonic 


l ahead of work change, see Fig. | 


2 behind work change, see Fig. | 


INTRODUCTION 


flow 


nozzle are 


ONE-DIMENSIONAI 


r I NHI 
pressible fluid through 


the relationships of heat addition and friction to a fluid in a chan 


ISENTROPIC equations of a com 


well known, as also are 


nel and the relationships across a shock wave. However, the re 

lationships for the flow of a compressible fluid between two areas 
in series with isentropic work addition or extraction have not 
been presented as such. This condition approximates numerous 
applications in the design and performance estimation of gas 


turbines and compressors 
DISCUSSION 


Fig. 1 is a diagram of the flow between two areas in series with 
an isentropic work addition or extraction between these two 
in terms of 


areas. For conservation of mass between two areas, 


the flow parameter WY 7/P4A, then 


WV T2/P2A2 = (WV 1T,/P1A1) V (T2/7 P,/P2) X 


and for an isentropic process, 
P\/P, = (1/1 


and since 


where a positive AJ is taken as power extraction, the continuity 


equation will become 








a 


ISENTROPIC 
WORK 
ADDITION 


OR 
yt) EXTRACTION 


A, A2 


Fic. 1. Schematic flow diagram 
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wT: [|(WVT 
PoAr | \ PoAe Jer 
wV 7, |(wV 7, wher 


PA, / PA, cr 
1 


(3) a ee 


Dy 


’ 


where (WY 72/P2A2)/(WYV To2/P2A2e),, is the ratio of the flow 
parameter at station 2 to the sonic flow at station 2 or is the ‘‘per 


cent of sonic flow”’ and & is the ratio of the per cent of sonic flow f 






op) 
= 
oO 
_j!. 
LL at station 2 to the per cent of sonic flow at station 1 
oO Fig. 2 shows the variation of © with area ratio (4;/A») and 0/on 
S change in specific output parameter (A7/7)). When the per The 
Oo cent of sonic flow is the same at both areas (e.g., both areas are “ 
Ww choked ), then ¥ is equal to 1 and the variation of the specific out ve 
ma put (A7/7,) is determined by the area ratio (.4;/A2) and js 
oO shown in Fig. 3. Values of = less than 1 correspond to the per a 
s cent of sonic flow larger at the first area than at the second area _— 
© (e.g., A; choked and A» unchoked). For = less than 1, as the _ 
a" specific output parameter increases |(A7/7)) increasing positive _ 
a for a given area ratio, the per cent of sonic flow at station 2 will st 
increase until it is the same at both areas (e.g., both areas are we 
choked). With a further increase in specific output, the per cent “e 
-.6 -4 -.2 ° a 4 6 of sonic flow at the second area will become larger than at the df 2 
TEMPERATURE DROP/inLET TEMPERATURE. 47/ first area (e.g., 42 remains choked and A; unchoked) and & will 
° bd become larger than 1. i 
Fic. 2. Sonic flow ratio variation with area ratio and specific Some applications of the equations of flow between two areas pert 
output in series with an isentropic work change between these areas time 
consist of the determination of the area ratio for maximum power ‘apy 
, output of a turbine, the matching and part load performance of W 
WY T2/P2A. = (WV 1)/P1A1) (A1/A2) X multistage turbines, the approximation of the equilibrium opera- eval 
{1 — (AT/T,)|°/2 — O/G—V) (1b) _ tion of a single spool turbojet engine (the turbine nozzle and tail- Her 
pipe nozzle represent the two areas) or of a twin spool engine or and 
From the isentropic channel flow equations, the flow parameter, other turbine and compressor combinations. nail 
WV T/PA, will be a maximum (critical value or “‘sonic flow’’) poin 
when the area is choked and is of tl 
: : -+ 
(WV T/PA),, = V 2y/[R(y + I) (2/¥ + DVA-YP (2) 
which is a constant for a given value of y. If the value of y is O 


the same at station 1 and station 2, then Evaluation of the Inertia Coefficients of the 
Cross Section of a Slender Body 

Arthur E. Bryson, Jr.* 

Harvard University, Cambridge, Mass laces 

















Pa January 25, 1954 

i = 

a 

© es coe ABSTRACT whe 
oa 

- - F T The determination of all the aerodynamic forces and moments (except and 
B= 2 drag) on a rigid slender body due to all possible motions of the body (except 

or T high frequency oscillations or step motions) reduces, for a certain large 

oO ae ee | class of slender configurations, to finding the six components of the sym whe 
= metric inertia coefficient tensor of the cross section of the body (see reference 1 
Ee '-+ 1). The evaluation of all six of these coefficients is given here in terms of the ; 
- coefficients of the Laurent Series expansion of the conformal mapping func un 
uj > we | tion that maps the cross section onto acircle. As examples the apparent ad- 

s ° ditional masses of the cross section of a wing-elliptic body-vertical tail con 

i mi | figuration and the apparent additional moment of inertia (which determines I 
ro) az i the roll damping) of the cross section of a configuration with » equal, equally fact 
x spaced fins are found 

oS “0 8 

uJ 

& = INTRODUCTION I 
= a a a a — THE FORCES on a two-dimensional body moving in any 

td a oe a eet | manner whatsoever through an infinite incompressible fluid, 

= ° 2 o 6 8 10 12 L4 where the fluid motion is irrotational and the fluid is at rest far 

- AREA RATIO, Ai / Az away from the body, reduces to the evaluation of the inertia co- 

one — a ae efficient tensor ; 
Fic. 3. Specific output variation for same per cent of sonic flow §——— 

out 


at both areas. * Assistant Professor, Division of Applied Science. 
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mix = Me; = —oF ¢; = ds j,k =1,2,3 (1 
on 
where 
density of the fluid 
g velocity potential due to unit velocity of the body 
in the x; direction (see for example Fig. 1 
g velocity potential due to unit velocity of the body 
in the x2 direction . 
¢ velocity potential due to unit angular velocity of - 
the body about the origin 
f = line integral around the contour of the body 
d = element of length along the contour 
)/0n = derivative normal to the contour 


The development of this fact is given in many places, for example, 
in Lamb.? 

One of the most interesting applications of this theory is to the 
so-called ‘‘slender body theory” (see, for example, reference 1) 
Direct use of the two-dimensional inertia coefficient concept in the 
slender body theory is apparently possible only for slender con 
figurations where the ‘‘fuselage’’ has two perpendicular planes of 
symmetry and all wings can be idealized as flat plates parallel 
to the axis formed by the intersection of the planes of symmetry 


of the fuselage. Thus reference 1 is too extravagant in its claims 


of generality. * 

The coefficients m1, M2, and m2 have the dimensions of mass 
per unit length and are often called ‘“‘apparent additional masses.” 
The coefficients 3 and m2; have dimensions of mass times length 
per unit length and the coefficient m3 has the dimensions of mass 
times length squared per unit length and can be thought of as an 
‘apparent additional moment of inertia.” 

Ward,*? Milne-Thompson,* and Summers® have shown how to 
evaluate m1, M2, and Me» in a simple manner by residue calculus. 
Here we wish to show a similar method for evaluating 3, m3, 
and m3 and we will repeat the method for evaluating the other 
coefficients simultaneously for completeness and a unified view- 
point. Two new examples will then be given showing the power 


of the method. 


I VALUATION OF 5; 


On the body surface 


dx» 1 
i (dx: | 
ds = —dx, : k= <2 (oO 
n 2 
0 ‘oe dx, — Xodx. 13 
hence 
My; + May = ipf oxdz (4) 
where z = x; + ix. and k Be 
and 
m33 = p/2 ft ¢3d(22) (5 
where 7 = the complex conjugate of z. 
The complex potential F; = ¢g, + 7~, where vy, is the stream 
function. Hence 
my + ima, = ip $ Fi.dz - if y;,dz) (6) 


Integrating the second integral of Eq. (6) by parts and using the 
fact that Oy; /O0s = Oy,/dn, we have 


F Vidz = — £ 2(d¢, On) ds 

are able to evaluate this integral as 
(-s ( 

{=< ; p <§ 

| —1Z)5 3 


Using Eq. (3) we 


* The author wishes to thank John R. Spreiter of the NACA for pointing 
out this fact to him 
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Fic. 1. The wing-elliptic body-vertical tail configuration 


where S = area of cross section and coordinate of center of 


area of cross section 
To find F; we map the cross section (in the z plane) conformally 


onto a circle of radius c in the ¢ plane, with no distortion at 


infinity. By Laurent’s theorem the mapping must then be of the 


form 
os . 
= Ns ha (8) 
n=O0¢ 
We then easily find that 
( c > ad, 
i - = 5 ~ / 5 = - T — n 
¢ ¢ n=O $ 
ic* sc — 
Ff = 241C — HC _ = —_ 3 (9) 
¢ ¢ ee 
" c? { ae 
F. = —iP.P ney ( ) =—1 3 : 
¢ f " ¢ 
where 
Am+n An a.\= 1 
h = >» = and : 10) 
an «4 co a = () form > 1 


where P.P. means “principal part of,” i.e., only the terms of the 





Laurent series expansion containing negative powers of ¢. From 
the residue theorem then 
(a, —< \1 
GF Fidz = 2rix —i(a; + c*); k= <2 (11 
Ff andy 13 
Substituting Eqs. (7) and (11) into Eq. (6) we have finally: 
c—a-— S$ 
my, + ime (c2 + a) — i(S/2e 
Mio + 1M2 | = 2rp : R (12) 
GU» UG» ‘ 
M3 + 1Mo3 l =: = Ga 
ra 2 9°. 
a ” 1 ¢ wi 


0, then my, = moe and my. = 0 


It is evident that if a; = 
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Fic. 2. The circular body with » equal, equally spaced fins 


The evaluation of ms; follows similar lines: 
p : 
M22 = | [fF d(zz) — if y d( 22 (13) 
but 


Oe: 
SL vid(2z) = — $23 eds 


on 
= (1/2) £'24d(=2) (14) 
= () 
Also, on the circle, ¢¢ = c?, so, on the body contour, 
ss = f(¢)f(c?/¢) (15) 
> * 


n=—o@ § 


where }, was defined in Eq. (10). It is easily shown from Eq. 


(10) that b_, = b,/c?", hence on the circle: 


by bag” 
22 = bb + a (*s cai e ) (16) 
n=1 $ 7 


if 


Substituting Eqs. (14) and (16) into Eq. (13) we have: 


. x 4 @ 
ip bn J Mbmo™ mbm | 
ms = — 9 ft jo +n 7  » a qm vine em +1 ( dj 
« n=1$ m=1 ¢c $ 
Nbnbn oe 
=p >, (17) 
-2n 
n=l ¢ 


Thus with Eqs. (12) and (17) we have all six components of the 
symmetric tensor ;, expressed in terms of the Laurent series 
coefficients of the conformal transformation of the cross section 
onto a circle through Eqs. (8) and (10). 


1—THE WING-ELLIPTICAL Bopy-VERTICAL TAIL 


CONFIGURATION 


EXAMPLE 


A slight generalization of the cross section considered in refer- 
ence 1 is the ellipse with two equal opposing fins and two unequal 
fins at right angles to the equal fins (see Fig. 1) (with & = a, 


SCIENCES JUNE, 1954 


this might be the cross section of a delta wing fighter with an elliy 
tical fuselage, highly swept delta wings and a vertical tail sur. 
face). 

The conformal mapping of this cross section onto a circle 


the ¢ plane is: 


a? — }? 
zs=w + 
472 
(a + 0)? |* h—f (h +f)? 
w + = a a) — & 18 
tr 2 F 16¢ 
where 
k = o6 — [(a + b)?/40 
h? = k? + [rm + (a + b)?/471]? 
f? = k? + [n+ (a b)?/4r 
and 
0 = 1/2s+ Vs? + a? — b 
mr = 1/24; Vt? — a2 + Bb? 
t2 = 1/2(t2 + Vt? — a? + BD? 


The first few terms of the Laurent series expansion are 


g=¢(+(h — f)/2 + 1/16¢ [(h + f)? — 8k? — 

Sh(a + b + (1/7 19 
and the radius of the circle in the ¢ plane isc = ( + f)/4. The 
area of the cross section is S = rab 


Therefore from Eq. (12), we have: 


My = m3 = 0 
my = mwp(k? + b? (20 
Moe = rp(4c? — k? — 2ab — Bb? 


The values of #2; and m3; are difficult to find in view of the fact 
the general term, a,, of the Laurent series, Eq. (8), for this case 
is a fairly complicated expression. The correct value for m, 
when a = 3, i.e., when the ellipse becomes a circle, is given in 
reference 6. 
EXAMPLE 2—CONFIGURATION WITH ” EQUAL, EQUALLY-SPACED 
FINS 
Miles’ has given the conformal mapping of a circle with n 


equal, equally spaced fins (see Fig. 2) onto a circle as 
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The apparent additional moment of inertia for a con- 
figuration with » equal, equally spaced fins 
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The Laurent expansion of Eq. (21) is 


s= f+ (2/n)[(c" — a*)/e— J] + O(1/F2"7! (22 
ence for n > 3, a, = O, and 
i = Ml = 71 = () 
bY" 2 + (a*/b 2 \4 ‘ 
1) = Mo = 2Zrp _ 23 
iy} } | 9 > { “v0 
Jready given in reference 7 
If a = 0 (no body) the Laurent expansion becomes simply a 
inomial series of the form 
—~ @ ril + (2/n i 
:= S ; where d,,, = — - (24 
ae ' rijl + (2/n) — m) m 
Substituting this expression for dm,—; into Eq. (10) we have 
fae, laieetas sk 
a 
— 
) 0 c 
ont - 9\")2 
= sin? rii- ri 2x 
r? n n 
> r{(] — (2/n) + m|T[—(2/n) + m ao 
~) 
Pigs r(/ + 1 + m)m! 
where use has been made of the relation 1/I'(1 — z) = (sin wz 


r) T(z 
This series is a hypergeometric series with argument unity 


which possesses a closed form sum,* so that Eq. (25) can be written 


—] _ 2a } P r{l — (2/n)]) , 
= sin riil4 = — (26 
7 n n {i + (2/n) + 1 


! (r[l — (2/n)|\? 


1+ (2/n)|? (r[l + (2 n)]§ 


99> 


(<4 


> 
l 


since c = b/2°/" 


1, 2, 4, 


ready found for these cases (reference 9), namely, 


Eq. (27) can be summed for ©, yielding values al 


97/128 


7/8 }. 28 
M33 = pb : n= 


ae r. 


W/m 


©, use is made of the asymptotic expression for 
becomes a Riemann 


where, for n — 
the I function, so that the series in Eq. (27 
Zeta Function: 


, 
‘one 
Considered as a function of the variable — = 1 + (8/m), this 
function has a simple pole with residue unity at — = 1 (see refer- 
ence & In fact 
8 1 8} 1 
*. = 1+ + 0 ( , asnu— 
hall n n? 
where > = Euler’s constant. Thus 
- 8 log 2 1 | 
Nieg = pb} 1 — + 0 = asn— (29 
2 n n* | 
The value for » = 3 was calculated to be m3; = 0.533 pb‘. 


Interpreted in terms of the roll damping (see reference 1) of 
a slender rocket or missile with n equal, equally spaced fins, we 


Nave 
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C; i 
(Cip)n (m 


A plot of this function is given in Fig. 3 
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Erratum —‘‘Analysis of Stiffened Curved Panels 
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[’ HAS BEEN CALLED to our attention that there is a misprint 


in Eq. (7) of our above-mentioned article. The exponent 1.125 
] 


should be 1.25 
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On Solving Subsonic Unsteady Flow Lifting 
Surface Problems by Separating Variables 


J. W. Miles 

Associate Professor of Engineering, University of California, Los 
Angeles 

February 8, 1954 


T° A RECENT PAPER dealing with the subsonic unsteady flow over 
a lifting surface, Professor Kuessner! has advocated the de 
velopment of ‘‘new orthogonal coordinate systems and new corre- 
sponding functions”’ in order to ‘‘put the problem on the strong 
shoulders of the mathematician.”’ It is possible that this state 
ment may lead to undue optimism, and it therefore seems perti 
nent to point out that the problem of solving the Helmholtz 


equation 
V*o + x*o@ = O | 


by separation of variables has been definitively studied by Robert 


son? and Ejisenhart,? who have shown that separation is pos- 


sible only in eleven (Euclidean) coordinate systems, viz., (1) rec- 


tangular, (2) circular cylinder, (3) elliptic cylinder, (4) parabolic 


evlinder, (5) spherical, (6) conical, (7) parabolic, (8) prolate 


spheroidal, (9) oblate spheroidal, (10) ellipsoidal, and (11) para 


boloidal. 
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Of these eleven systems only the three discussed by Kuessner 
define lifting surfaces of practical interest, viz., elliptic cylinder 
(infinite, flat ribbon; i.e., two dimensional wing), oblate sphe- 
roidal (circular planform), and ellipsoidal (elliptic planform) 
To these might be added the coordinates of the parabolic cylin- 
der, which may be used in the analysis of a semi-infinite plane, 
but such a lifting surface is of rather limited, practical interest in 
subsonic flow. Accordingly, it must be inferred that any expec- 
tation of applying the method of separation of variables to the 
solution of new, lifting surface problems in unsteady, subsonic, 
compressible flow is unwarranted 

It should be added that Laplace’s equation [x = 0 in (1), corre- 
sponding to incompressible flow] is separable in additional co- 
ordinate systems such as: the infinity of two dimensional systems 
obtainable by conformal transformation; additional three di- 
mensional systems following from their two dimensional counter- 
parts by rotation;* and confocal cyclides® (e.g., toroidal coordi- 
nates); however, it appears that none of these lead to new (wing) 
surfaces of practical interest 

The foregoing remarks are limited to Euclidean space, and 
other systems of coordinates may be developed for non-Euclidean 
spaces. Thus, in supersonic flow a modified Lorentz transfor- 
mation of the equations of motion to a coordinate system mov- 
ing with the lifting surface leads to a hyperbolic space having the 
differential metric (ds)? = (dx)? — (dy)? — (dz)?. Unsteady 
flow problems in this space then lead to the ‘hyperbolic Helm- 
holtz’’ equation 


dir — hyy — bez + Kd = 0 (2) 


This equation also is separable in eleven coordinate systems, 


seven of which differ from 1-11 above. Of these seven, Robin- 
son’s hyperboloido-conal coordinates’ (having the delta wing as 


one surface) are perhaps the most useful in practical application 
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On Down Flapping of a Stalled Helicopter 
Blade* 


W. Z. Stepniewski 

Director of Aerodynamics, 
Morton, Pa. 

February 14, 1954 


Piasecki Helicopter Corporation, 


—— NOTE is intended to indicate the possibility of an excessive 
down-flapping motion which may be developed by a stalled 
helicopter blade. Since such motion may occur either in flight 
maneuvers or at landing, its practical significance may be of 
interest to those engaged in design and operation of helicopters. 

In order to give some feeling as to the possible order of magni- 
tude of the down-flapping of stalled blade, and at the same time 

* Opinions expressed in this note are those of the writer alone and not 
necessarily those of Piasecki Helicopter Corporation. 
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outline only the most important physical parameters influencing 
it, the following simplifying assumptions are made: 

(1) An articulated blade of constant chord c, having it flap- 
ping hinge at the rotor axis and belonging to a helicopter moving 
at a moderatef tip speed ratio u, becomes completely stalled jy 


the up-wind position, i.e., at an azimuth angle yo = 7 
(2) At this moment of stalling, the blade reaches its maximum 
up-flapping angle 8) = 8,, which in turn, requires that 8, = 0 
(3) The whole blade which became stalled at Yo = z, remains 


in this condition at least up to Y = 27, and throughout the whok 
azimuth angle  < y < 27m the lift coefficient of the blade ele. 
ments remains constant being c,) = c;, = constant; where ¢, js 
a representative value of c of the ‘‘horizontal’’ portion of the 
c, = f(@) curve in the region of a fully developed stall 
Knowing boundary conditions as given by yo, Bo and 8», and 
accepting the above assumptions, it is rather easy to find an 
equation describing the blade flapping motion from the assumed 
moment of stall (yo = w) until it reaches its down-wind position 


(y = 2m). The greatest simplification of the above task stems 
from the assumption of c; = c;, = constant, as it makes the 


aerodynamic moment (./4) about the flapping hinge independ 
ent of the blade flapping velocity, which in turn, permits one to 
obtain a very simple expression for .\/ 4, easily derived from the 
elementary aerodynamic moment; 

dM, = (1/2)pcR%c,,R%w2(x + uw sin ~)? x dx (1 
where p is air density (in slugs/cu.ft.) ¢ is blade chord (in ft 
R is blade radius (also in ft.), w is angular velocity (rad./sec 
of the rotor, x is the nondimensional blade station, and all other 
symbols as previously defined 

Neglecting the reversed flow 
ment experienced by the blade as a whole, can easily be obtained 
(1) from zero to the effective station x, 


region, the aerodynamic mo 


by integrating Eq 


» 


Ma = (1/2)pcRw*c;,[(1/4)xe4 + (2/3)x.34 sin Y + 
(1/2)xe%4? sin? y] (2 


or, assuming for instance xp = 0.97, Eq. (2) becomes 


Ma = (1/2)pcR4w*c;,(0.221 + 0.6074 sin y 4 


0.470p2 sin? y¥ (2a 

Neglecting the influence of gravity, the basic equation of the 
blade flapping motion can be written as follows: 

B+ w8 = Ma4/I; (3 

where 7; is the blade moment of inertia about the flapping hinge 


Introducing Eq. (2a) into Eq. (3), the following is obtained: 


B+ ws = [(pcR*c;, ) 277] w?(0.221 + 0.607 sin Y + 
0.470? sin? y (4 


It should be noticed that on the right-hand side of Eq. (4), 
the factor pcR'c;,/2I; resembles closely, the Lock number, y = 
pcR'a/2I;,, with the only difference that the lift curve slope, a, 
is replaced by the lift coefficient of the stalled blade « To 
this new factor the symbol ys will be assigned as normally used 


for Lock number: 


¥s = [(pcR4c1, ) 2T 7) (oO) 


It has been assumed that the blade becomes stalled in its up- 
wind position, and its motion is investigated starting from this 
moment. This means that zero of the time scale corresponds to 
y = 7m, and the azimuth angle during the investigated period of 


blade motion can be expressed as: 
y=7+ a (6) 


Introducing Eqs. (5) and (6) into Eq. (4), the following is ob- 


tained: 


B+ wB = yew?(0.221 — 0.607u sin wt + 0.470u? sin? wt) (7 


t Tip speed ratios up to about 4 = 0.35 may be considered as moderate 


as they permit to neglect the influence of reversed velocities 
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In order to solve Eq. (7), it is first operated with Laplace trans- 


formations, becoming: 


23 — sa(0O) — B(O) + ws = 
().221 0.607 pw 0.235y? 0.2352 
yw? + — (S 
s s2 + w? , s2 + (Qu)? 
Taking into consideration boundary conditions at ¢ = 0, 
3) = B, and B = 6, = 0, Eq. (8) can be rewritten as follows: 
§ = Bp, = i yee 
52 + wy? 
40.221 + 0.235p? 0.607 pw 
— T 
/ s(s? + w?) (s? + w*)(s* + w*) 
0 2352s { 


(Sa ) 


[s* + (2w)*](s? + w?)4 


All terms on the right side of Eq. (Sa) are such that their inverse 
transformation can easily be performed with standard tables of 


Laplace transforms (see, for instance, Appendix III of reference 


1), and the desired expression for 8 is readily obtained, 
8 = B, cos wt + ys[(0.221 + 0.235u?)(1 — cos wt) — 
0.303u(sin wt — wt cos wt) — O.078u2(cos wt — cos 2wt)| (9 


Substituting for w/ its value from Eq. (6), Eq. (9) can be re 


written in terms of the azimuth angle: 


B) ‘ 9 7 
3 = —B, cos yy + ¥s}0.221(1 + cos y) + 


0.303u|sin Yy — (~¥ — mw) cos ¥] 4 

p?(0.235(1 + cos Y) + O O78(cos y + cos 2y)]; (9a 

In the down-wind position (y = 27), which is of special interest 
from the practical point of view, the flapping angle becomes: 

Bor = —Ba + y;(0.442 — 0.950 + 0.6267) (10) 


An inspection of Eq. (10) will clearly indicate that the amount of 
down-flapping of stalled blade in the aft position depends, first 
f all, on the flapping angle which it had in the up-wind position 
This means that of control displacements, 
rotor loading, r.p.m., ete., which can result in a high up-wind 
flapping angle of the blade at the moment of stalling, will directly 
contribute to the low flapping of the stalled blade in its down 


any combination 


wind position 

The second term on the right side of Eq. (10) remains positive 
for all moderate values of yw, but its magnitude decreases with 
increasing tip speed ratios. The character of this second term 
means that, in general, it tends to decrease the amount of down- 


flapping in the down-wind position of the blade. Since the 
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value of the whole second term is proportional to yx, it is clear 
that of ys (the higher blade 


moment factors remain the 


the lower is this value flapping 


of inertia when all other same) the 


smaller restraining influence this term will have on the down 
flapping of the stalled blade 

In order to develop some feeling as to the order of magnitude of 
the ys values which can be encountered in practice, one should 
remember that y, differs from Lock number by replacing the lift 
slope a in this latter expression by ¢ 

Considering that c, is of the order of magnitude of 0.9 or 1.0, 
while slope of the lift coefficient is usually assumed as 
that 


Lock number 


= » 


to about one 


sixth of the Lock 
equal to 4 or 5 may be considered as fairly representative of actual 


it becomes clear the value of ys, will amount 


Since in addition a number 
construction blades, y, = 0.7 may be used as a practical value 


Fig. 1 gives 6, = {(8,) for the above derived value of y; = 
0.7. In addition, some feeling as to the influence of tip speed 
ratio can be obtained from this figure as the amount of down 
flapping is shown for three values of u 

A glance at Fig. | will confirm the previously mentioned state 
ment that any combination of conditions which may lead to stall 
ing of the blade at high flapping angle in its up-wind position, 
may produce an appreciable amount of negative flapping when 


the blade reaches the down-wind position 
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A Paradox in the Theory of Impact 


C. Riparbelli 

Associate Professor, Graduate School of Aeronautical Engineering, 
Cornell University, Ithaca, N 

February 19, 1954 


A 


strain in % 
end impact—i.e., at whose free end a constant velocity is imposed 


PARADOX appears in the theory of longitudinal impact 
More precisely, on the subject of the distribution of plastic 
i semi-infinite bar of uniform cross section undergoing 
for a certain duration. Most of experimental data are about 
soft copper, so our discussion will be limited to that material 

Here are the terms of the paradox 

(a) Assuming the static stress strain diagram to be 
neglecting the time-rate effect on the 


valid in 
dynamic conditions (1.e, 
stress-strain relation), a distribution of plastic strain along the 
bar, function of time, is predicted.! 
the experimental data on such distri 


It is seen that it follows 


at least 
bution, both on measure of residual strain after impact? 


qualitatively 
and by 
direct measure of strain during the impact 

(b) A group of data are available which do not fit the hypoth 
Namely 
stress-strain 


esis at (a) 
(1) The 


Direct measure of the stress in impact conditions indicates ulti 


relation is highly time dependent 


mate stresses up to three times the static one 

(2) The velocity of propagation of a strain-raising pulse is 
always the one correspondent to elastic behavior of the material, 
even if the bar is prestressed high in the plastic range*: '® and even 
if it is undergoing plastic flow while the pulse travels.'! This 
velocity is measured at the front of the pulse 

(3) The traveling pulse, defined at the beginning in its length 
by the impact duration, does not remain of constant length while 


traveling, but it spreads, while its maximum amplitude de 
creases. !® 1 
(4) The stress measured directly at the end under impact (by 


the stress 


measure of the hammer’s acceleration) shows that 
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The 


plot of the stress and strain pairs on the stress-strain plane gives 


raise is immediate, while the plastic strain follows in time.’ 


points much higher than those corresponding to the static rela 
tion.’ 

The apparent contradiction between (a) and (b) is eliminated 
by introducing the plastic flow,’ i.e., the time rate of plastic 
deformation as an essential parameter. The plastic strain re 
sults correctly as the resultant of plastic flow during the history 
of the impact.'?) '* The time rate of plastic flow can be com 
puted with sufficient approximation according to a relation of 
quasi-viscosity within a limited range of time rates of strain, 
above which the stress increases with the time rate of strain but 
less than linearly A finite stress corresponds to time rates of 
strain tending to infinity.” § 

At every point of the bar the first stress raise is elastic. The 
plastic strain grows in time while the elastic strain decreases ac 
cording to the end conditions. The plastic flow is stress-reliev 
ing, so that the region of the bar nearer to the end under impact 
after sufficient duration of the impact is found to undergo a stress 
lower than the static stress corresponding to the strain present 
(i.e., according to the strain hardening). This region of the bar 
moves with uniform velocity and uniform stress, as predicted in 
the theory according to (a).'* 

It is also seen that within the region of quasi-viscous relation 
between the stress and the time rate of strain, the lines of equal 
strain rate are equal to the static stress-strain diagram, but dis 
placed in the direction of Hooke’s line.*»7 | Therefore the deriva- 
tive of such diagrams, i.e., ratio of stress increment to increment 
of plastic strain, versus plastic strain alone, is independent of the 
stress present. It is concluded that within the limit of the linear 
relation between the stress and the time rate of strain, it is im 
possible to conclude anything about the stress by measuring only 
the distribution of plastic strain 
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espa interest has been manifested in the problem of the 
effect of leading-edge bluntness on the pressures on a flat 
high Mach 
affects the pressures.'! Data bearing on this subject have been 
Mach 
In this investigation, pressure distributions 


plate at Numbers when the boundary layer also 


obtained at the Langley 1l-in. hypersonic tunnel at a 
Number of 6.86 
were obtained along a 20° wedge probe with various leading-edge 
The thicknesses 
were varied by making a plane cut normal to the surface on which 


thicknesses at a constant free stream pressure 
the pressures were measured. These data are presented in Fig 


1 as the nondimensional pressure rise (Ps — P;)/P, against the 
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ratio of length to leading-edge thickness, where P» is the local 
static pressure and 


an inviscid rotational flow solution computed 


the pressure in the undisturbed stream 
For comparison 
gr 
4 wedge leading edge is included 


blunt 


iphically, in a manner similar to that given by Meyer,? for a 
t plate with a 


flat The 
vedge leading edge of the theoretical model is of such an angle 
isto give Mach Number 1 behind the shock (about 43° for J, = 
6.86 Chis curve, of course, depends only on the dimensionless 
rameter x/f. From the experimental data, it is readily appar- 
nt that, unlike the flow 
Reynolds Number forms a separate curve though all the curves 


inviscid solutions each leading-edge 
ire of similar shape 

From the faired results in Fig. 1, Fig 
which, in effect, presents the pressure rise as a function of the 


2 has been prepared 


leading-edge thickness where R,, is the Reynolds number based 
on conditions in the undisturbed stream with leading-edge thick 
ness as the characteristic length. The cross-plotted experimental 
lata have been faired to zero thickness to give the effect of the 
boundary laver on the surface pressures unmodified by leading 
edge pressure effects due to thickness. The values obtained for 
' = 0 are in good agreement with pressures obtained with the 
issumption that the ordinary displacement thickness of the bound 
iy layer is the viscous layer which displaces the free stream 
flow and thereby alters the pressures over a plate As an as 
sumed first approximation to blunt leading-edge effects where 
boundary-layer effects are already present, the inviscid rotational 
flow solution shown on Fig. | has been plotted on Fig. 2 increased 
by the pressure rise indicated by the extrapolated experimental 
data at Re = O. At Rez = 30,000 and Re, = 2,000 the experi 
mental data indicate an increase in pressure of about 100 per 
cent over the zero thickness case while the rotational flow curve 
for these conditions indicates an increase of about 120 per cent 
At this stage, this is reasonable agreement considering the scatter 
of the experimental data and the simple assumption involved 
Fig. 3 presents experimental pressure data from various flat 
plates in addition to the 20° wedge probe in comparison to results 
from a displacement method® including the effects of the self 
induced pressure gradient on the variation of 6* with x variable 
dM/dx). The 
component of velocity and the pressure gradient in the vertical 
direction negligible V R, is that 
for which the theory developed by Lees‘ predicts a linear rela 
tionship with pp». 


> 


sented in Fig. 3 have been presented previously.* 


method assumes an insulated plate with the 


The abscissa on this figure 1 


Most of the experimental pressure data pre- 
5 The dis 
placement method after a solution by iteration starting from the 
result of a parabolic 6* = 6*(x) variation (d. /dx = 0) gives re 
pi(1/V R,.). While 


sults which are essentially linear in p. = 
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the result of a parabolic displacement thickness variation does 
not give a linear result when plotted in this fashion, it does not 
accurate curve for the range of 
Reynolds number shown on Fig. 3 at this Mach Number Also 


shown on Fig. 3 is the result given by the displacement method 


differ greatly from the more 


») 


when corrected for leading-edge thickness according to Fig. 2 
These results indicate that, with the size of models currently 


being tested in hypersonic tunnels, quite large pressure effects 
due to leading-edge thickness are possible even with models 
whose leading edges are generally considered quite sharp. Such 


into account in the experimental deter 


effects must be taken 
mination of the effects of the boundary layer on the pressures on 


a flat plate 
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